Class PCGSolver

java.lang.Object
ffx.potential.nonbonded.pme.PCGSolver

public class PCGSolver extends Object
Parallel pre-conditioned conjugate gradient solver for the self-consistent field.

This solves the linear system A x = b using an iterative approach where

A: is an N x N symmetric, positive definite matrix

B: is a known vector of length N.

x: is the unknown vector of length N.

For the AMOEBA SCF, the linear system Ax = b is usually denoted: C u = E_direct
where C = [alpha^-1 - T]
u are the induced dipoles.
E_direct is the direct field from permanent multipoles.

The matrix alpha^-1 is the inverse of the N x N diagonal polarizability matrix. The matrix T is the N x N matrix that produces the field due to induced dipoles.

Initialization:
1) Compute the residual where x_0 is either u_direct or u_mutual (a pcg guess):
r_0 = b - A x_0
r_0 = E_direct - C u_0
r_0 = E_direct - [alpha^-1 - T] u_0
r_0 = (u_direct - u_0) * alpha^-1 - E_u_0
2) Compute the preconditioner: z_0 = M^1 r_0
3) Compute the conjugate: p_0 = z_0
4) Initial loop index: k = 0

Then loop over:
1) Update the step size: alpha = r_k dot z_k / (p_k A p_k)
2) Update the induced dipoles: u_k+1 = u_k + alpha * p_k
3) Update the residual: r_k+1 = r_k - alpha * A p_k
4) Check for convergence of r_k+1.
5) Update the preconditioner: z_k+1 = M^-1 r_k+1
6) Update the step size: beta = r_k+1 dot z_k+1 / (r_k dot z_k)
7) Update the conjugate: p_k+1 = z_k+1 + beta * p_k
8) Update the loop index: k = k + 1

Since:
1.0
Author:
Michael J. Schnieders
  • Field Details

    • DEFAULT_CG_PRECONDITIONER_CUTOFF

      public static final double DEFAULT_CG_PRECONDITIONER_CUTOFF
      A preconditioner cut-off of 3 to 4 Angstroms generally works well.
      z = M^(-1) r = polarizability * (E_r + scale * r)
      The cutoff is applied when computing E_r (the field due to polarizability * residual).
      See Also:
    • DEFAULT_CG_PRECONDITIONER_EWALD

      public static final double DEFAULT_CG_PRECONDITIONER_EWALD
      An Ewald coefficient of 0 indicates use of Coulomb's law with no damping.
      z = M^(-1) r = polarizability * (E_r + scale * r)
      Ewald or Coulomb is used when computing E_r (the field due to polarizability * residual).
      See Also:
    • DEFAULT_CG_PRECONDITIONER_SCALE

      public static final double DEFAULT_CG_PRECONDITIONER_SCALE
      The scale factor is applied to the diagonal terms of the preconditioner.
      z = M^(-1) r = polar * (E_r + scale * r)
      The "scale * polar * r" term is the diagonal part of M^(-1)
      See Also:
  • Constructor Details

    • PCGSolver

      public PCGSolver(int maxThreads, double poleps, ForceField forceField, int nAtoms)
      Constructor the PCG solver.
      Parameters:
      maxThreads - Number of threads.
      poleps - Convergence criteria (RMS Debye).
      forceField - Force field in use.
      nAtoms - Initial number of atoms.
  • Method Details

    • allocateLists

      public void allocateLists(int nSymm, int nAtoms)
      Allocate storage for pre-conditioner neighbor list.
      Parameters:
      nSymm - Number of symmetry operators.
      nAtoms - Number of atoms.
    • allocateVectors

      public void allocateVectors(int nAtoms)
      Allocate PCG vectors.
      Parameters:
      nAtoms - The number of atoms.
    • getPreconditionerCutoff

      public double getPreconditionerCutoff()
      Get the preconditioner cutoff.
      Returns:
      The cutoff in Angstroms.
    • getPreconditionerEwald

      public double getPreconditionerEwald()
      Get the preconditioner Ewald coefficient.
      Returns:
      The Ewald coefficient.
    • getPreconditionerScale

      public double getPreconditionerScale()
      Get the preconditioner matrix diagonal scale factor.
      Returns:
      The unitless scale factor.
    • getPreconditionerMode

      public String getPreconditionerMode()
      Get the preconditioner mode.
      Returns:
      The mode.
    • getPreconditionerLists

      public int[][][] getPreconditionerLists()
      Neighbor lists when applying the preconditioner.
      Returns:
      Storage for the preconditioner neighbor lists.
    • getPreconditionerCounts

      public int[][] getPreconditionerCounts()
      Number of neighbors when applying the preconditioner.
      Returns:
      Storage for the preconditioner counts.
    • init

      public void init(Atom[] atoms, double[][][] coordinates, double[] polarizability, double[] ipdamp, double[] thole, boolean[] use, Crystal crystal, double[][][] inducedDipole, double[][][] inducedDipoleCR, double[][] directDipole, double[][] directDipoleCR, AtomicDoubleArray3D field, AtomicDoubleArray3D fieldCR, EwaldParameters ewaldParameters, double dieletric, ParallelTeam parallelTeam, IntegerSchedule realSpaceSchedule, PMETimings pmeTimings)
    • scfByPCG

      public int scfByPCG(boolean print, long startTime, ParticleMeshEwald particleMeshEwald)