Class PCGSolver
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 Summary
Modifier and TypeFieldDescriptionstatic final double
A preconditioner cut-off of 3 to 4 Angstroms generally works well.static final double
An Ewald coefficient of 0 indicates use of Coulomb's law with no damping.static final double
The scale factor is applied to the diagonal terms of the preconditioner. -
Constructor Summary
ConstructorDescriptionPCGSolver
(int maxThreads, double poleps, ForceField forceField, int nAtoms) Constructor the PCG solver. -
Method Summary
Modifier and TypeMethodDescriptionvoid
allocateLists
(int nSymm, int nAtoms) Allocate storage for pre-conditioner neighbor list.void
allocateVectors
(int nAtoms) Allocate PCG vectors.int[][]
Number of neighbors when applying the preconditioner.double
Get the preconditioner cutoff.double
Get the preconditioner Ewald coefficient.int[][][]
Neighbor lists when applying the preconditioner.Get the preconditioner mode.double
Get the preconditioner matrix diagonal scale factor.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) int
scfByPCG
(boolean print, long startTime, ParticleMeshEwald particleMeshEwald)
-
Field Details
-
DEFAULT_CG_PRECONDITIONER_CUTOFF
public static final double DEFAULT_CG_PRECONDITIONER_CUTOFFA 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_EWALDAn 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_SCALEThe 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
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
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
-