List of usage examples for org.apache.commons.math3.optimization SimpleVectorValueChecker SimpleVectorValueChecker
public SimpleVectorValueChecker(final double relativeThreshold, final double absoluteThreshold)
From source file:ffx.potential.nonbonded.ParticleMeshEwaldCart.java
/** * ParticleMeshEwald constructor.//w ww .j a va2 s. c om * * @param atoms the Atom array to do electrostatics on. * @param molecule the molecule number for each atom. * @param forceField the ForceField the defines the electrostatics * parameters. * @param crystal The boundary conditions. * @param elecForm The electrostatics functional form. * @param neighborList The NeighborList for both van der Waals and PME. * @param parallelTeam A ParallelTeam that delegates parallelization. */ public ParticleMeshEwaldCart(Atom atoms[], int molecule[], ForceField forceField, Crystal crystal, NeighborList neighborList, ELEC_FORM elecForm, ParallelTeam parallelTeam) { this.atoms = atoms; this.molecule = molecule; this.forceField = forceField; this.crystal = crystal; this.parallelTeam = parallelTeam; this.neighborList = neighborList; this.elecForm = elecForm; neighborLists = neighborList.getNeighborList(); permanentSchedule = neighborList.getPairwiseSchedule(); nAtoms = atoms.length; nSymm = crystal.spaceGroup.getNumberOfSymOps(); maxThreads = parallelTeam.getThreadCount(); polsor = forceField.getDouble(ForceFieldDouble.POLAR_SOR, 0.70); poleps = forceField.getDouble(ForceFieldDouble.POLAR_EPS, 1e-5); if (elecForm == ELEC_FORM.PAM) { m12scale = forceField.getDouble(ForceFieldDouble.MPOLE_12_SCALE, 0.0); m13scale = forceField.getDouble(ForceFieldDouble.MPOLE_13_SCALE, 0.0); m14scale = forceField.getDouble(ForceFieldDouble.MPOLE_14_SCALE, 0.4); m15scale = forceField.getDouble(ForceFieldDouble.MPOLE_15_SCALE, 0.8); } else { double mpole14 = 0.5; String name = forceField.toString().toUpperCase(); if (name.contains("AMBER")) { mpole14 = 1.0 / 1.2; } m12scale = forceField.getDouble(ForceFieldDouble.MPOLE_12_SCALE, 0.0); m13scale = forceField.getDouble(ForceFieldDouble.MPOLE_13_SCALE, 0.0); m14scale = forceField.getDouble(ForceFieldDouble.MPOLE_14_SCALE, mpole14); m15scale = forceField.getDouble(ForceFieldDouble.MPOLE_15_SCALE, 1.0); } d11scale = forceField.getDouble(ForceFieldDouble.DIRECT_11_SCALE, 0.0); p12scale = forceField.getDouble(ForceFieldDouble.POLAR_12_SCALE, 0.0); p13scale = forceField.getDouble(ForceFieldDouble.POLAR_13_SCALE, 0.0); useCharges = forceField.getBoolean(ForceFieldBoolean.USE_CHARGES, true); useDipoles = forceField.getBoolean(ForceFieldBoolean.USE_DIPOLES, true); useQuadrupoles = forceField.getBoolean(ForceFieldBoolean.USE_QUADRUPOLES, true); rotateMultipoles = forceField.getBoolean(ForceFieldBoolean.ROTATE_MULTIPOLES, true); lambdaTerm = forceField.getBoolean(ForceFieldBoolean.LAMBDATERM, false); if (!crystal.aperiodic()) { off = forceField.getDouble(ForceFieldDouble.EWALD_CUTOFF, 7.0); } else { off = forceField.getDouble(ForceFieldDouble.EWALD_CUTOFF, 1000.0); } double ewaldPrecision = forceField.getDouble(ForceFieldDouble.EWALD_PRECISION, 1.0e-8); aewald = forceField.getDouble(ForceFieldDouble.EWALD_ALPHA, ewaldCoefficient(off, ewaldPrecision)); setEwaldParameters(off, aewald); reciprocalSpaceTerm = forceField.getBoolean(ForceFieldBoolean.RECIPTERM, true); String predictor = forceField.getString(ForceFieldString.SCF_PREDICTOR, "NONE"); try { predictor = predictor.replaceAll("-", "_").toUpperCase(); scfPredictor = SCFPredictor.valueOf(predictor); } catch (Exception e) { scfPredictor = SCFPredictor.NONE; } if (scfPredictor != SCFPredictor.NONE) { predictorCount = 0; int defaultOrder = 6; predictorOrder = forceField.getInteger(ForceFieldInteger.SCF_PREDICTOR_ORDER, defaultOrder); if (scfPredictor == SCFPredictor.LS) { leastSquaresPredictor = new LeastSquaresPredictor(); double eps = 1.0e-4; leastSquaresOptimizer = new LevenbergMarquardtOptimizer(new SimpleVectorValueChecker(eps, eps)); } else if (scfPredictor == SCFPredictor.ASPC) { predictorOrder = 6; } predictorStartIndex = 0; } String algorithm = forceField.getString(ForceFieldString.SCF_ALGORITHM, "CG"); try { algorithm = algorithm.replaceAll("-", "_").toUpperCase(); scfAlgorithm = SCFAlgorithm.valueOf(algorithm); } catch (Exception e) { scfAlgorithm = SCFAlgorithm.CG; } /** * The size of the preconditioner neighbor list depends on the size of * the preconditioner cutoff. */ if (scfAlgorithm == SCFAlgorithm.CG) { inducedDipolePreconditionerRegion = new InducedDipolePreconditionerRegion(maxThreads); pcgRegion = new PCGRegion(maxThreads); pcgInitRegion1 = new PCGInitRegion1(maxThreads); pcgInitRegion2 = new PCGInitRegion2(maxThreads); pcgIterRegion1 = new PCGIterRegion1(maxThreads); pcgIterRegion2 = new PCGIterRegion2(maxThreads); boolean preconditioner = forceField.getBoolean(ForceFieldBoolean.USE_SCF_PRECONDITIONER, true); if (preconditioner) { preconditionerCutoff = forceField.getDouble(ForceFieldDouble.CG_PRECONDITIONER_CUTOFF, 4.5); preconditionerEwald = forceField.getDouble(ForceFieldDouble.CG_PRECONDITIONER_EWALD, 0.0); } else { preconditionerCutoff = 0.0; } } else { preconditionerCutoff = 0.0; inducedDipolePreconditionerRegion = null; pcgRegion = null; pcgInitRegion1 = null; pcgInitRegion2 = null; pcgIterRegion1 = null; pcgIterRegion2 = null; } if (lambdaTerm) { /** * Values of PERMANENT_LAMBDA_ALPHA below 2 can lead to unstable * trajectories. */ permLambdaAlpha = forceField.getDouble(ForceFieldDouble.PERMANENT_LAMBDA_ALPHA, 2.0); if (permLambdaAlpha < 0.0 || permLambdaAlpha > 3.0) { permLambdaAlpha = 2.0; } /** * A PERMANENT_LAMBDA_EXPONENT of 1 gives linear charging of the * permanent electrostatics, which is most efficient. A quadratic * schedule (PERMANENT_LAMBDA_EXPONENT) also works, but the dU/dL * forces near lambda=1 are may be larger by a factor of 2. */ permLambdaExponent = forceField.getDouble(ForceFieldDouble.PERMANENT_LAMBDA_EXPONENT, 1.0); if (permLambdaExponent < 1.0) { permLambdaExponent = 1.0; } /** * A POLARIZATION_LAMBDA_EXPONENT of 2 gives a non-zero d2U/dL2 at * the beginning of the polarization schedule. Choosing a power of 3 * or greater ensures a smooth dU/dL and d2U/dL2 over the schedule. */ polLambdaExponent = forceField.getDouble(ForceFieldDouble.POLARIZATION_LAMBDA_EXPONENT, 3.0); if (polLambdaExponent < 0.0) { polLambdaExponent = 0.0; } /** * The POLARIZATION_LAMBDA_START defines the point in the lambda * schedule when the condensed phase polarization of the ligand * begins to be turned on. If the condensed phase polarization is * considered near lambda=0, then SCF convergence is slow, even with * Thole damping. In addition, 2 (instead of 1) condensed phase SCF * calculations are necessary from the beginning of the window to * lambda=1. */ polLambdaStart = forceField.getDouble(ForceFieldDouble.POLARIZATION_LAMBDA_START, 0.75); if (polLambdaStart < 0.0 || polLambdaStart > 0.9) { polLambdaStart = 0.75; } /** * The POLARIZATION_LAMBDA_END defines the point in the lambda * schedule when the condensed phase polarization of ligand has been * completely turned on. Values other than 1.0 have not been tested. */ polLambdaEnd = forceField.getDouble(ForceFieldDouble.POLARIZATION_LAMBDA_END, 1.0); if (polLambdaEnd < polLambdaStart || polLambdaEnd > 1.0 || polLambdaEnd - polLambdaStart < 0.3) { polLambdaEnd = 1.0; } /** * The LAMBDA_VAPOR_ELEC defines if intramolecular electrostatics of * the ligand in vapor will be considered. */ doLigandVaporElec = forceField.getBoolean(ForceFieldBoolean.LIGAND_VAPOR_ELEC, true); doNoLigandCondensedSCF = forceField.getBoolean(ForceFieldBoolean.NO_LIGAND_CONDENSED_SCF, true); /** * Flag to indicate application of an intermolecular softcore * potential. */ intermolecularSoftcore = forceField.getBoolean(ForceField.ForceFieldBoolean.INTERMOLECULAR_SOFTCORE, false); intramolecularSoftcore = forceField.getBoolean(ForceField.ForceFieldBoolean.INTRAMOLECULAR_SOFTCORE, false); } String polar = forceField.getString(ForceFieldString.POLARIZATION, "MUTUAL"); if (elecForm == ELEC_FORM.FIXED_CHARGE) { polar = "NONE"; } boolean polarizationTerm = forceField.getBoolean(ForceFieldBoolean.POLARIZETERM, true); if (polarizationTerm == false || polar.equalsIgnoreCase("NONE")) { polarization = Polarization.NONE; } else if (polar.equalsIgnoreCase("DIRECT")) { polarization = Polarization.DIRECT; } else { polarization = Polarization.MUTUAL; } String temp = forceField.getString(ForceField.ForceFieldString.FFT_METHOD, "PJ"); FFTMethod method; try { method = ReciprocalSpace.FFTMethod.valueOf(temp.toUpperCase().trim()); } catch (Exception e) { method = ReciprocalSpace.FFTMethod.PJ; } gpuFFT = method != FFTMethod.PJ; if (lambdaTerm) { shareddEdLambda = new SharedDouble(); sharedd2EdLambda2 = new SharedDouble(); } else { shareddEdLambda = null; sharedd2EdLambda2 = null; lambdaGrad = null; lambdaTorque = null; vaporCrystal = null; vaporLists = null; vaporPermanentSchedule = null; vaporEwaldSchedule = null; vacuumRanges = null; } if (logger.isLoggable(Level.INFO)) { StringBuilder sb = new StringBuilder("\n Electrostatics\n"); sb.append(format(" Polarization: %8s\n", polarization.toString())); if (polarization == Polarization.MUTUAL) { sb.append(format(" SCF Convergence Criteria: %8.3e\n", poleps)); sb.append(format(" SCF Predictor: %8s\n", scfPredictor)); sb.append(format(" SCF Algorithm: %8s\n", scfAlgorithm)); if (scfAlgorithm == SCFAlgorithm.SOR) { sb.append(format(" SOR Parameter: %8.3f\n", polsor)); } else { sb.append(format(" CG Preconditioner Cut-Off: %8.3f\n", preconditionerCutoff)); sb.append(format(" CG Preconditioner Ewald Coefficient:%8.3f\n", preconditionerEwald)); } } if (aewald > 0.0) { sb.append(" Particle-mesh Ewald\n"); sb.append(format(" Ewald Coefficient: %8.3f\n", aewald)); sb.append(format(" Particle Cut-Off: %8.3f (A)", off)); } else { sb.append(format(" Electrostatics Cut-Off: %8.3f (A)\n", off)); } logger.info(sb.toString()); } if (gpuFFT) { sectionThreads = 2; realSpaceThreads = parallelTeam.getThreadCount(); reciprocalThreads = 1; sectionTeam = new ParallelTeam(sectionThreads); realSpaceTeam = parallelTeam; fftTeam = new ParallelTeam(reciprocalThreads); } else { boolean concurrent; int realThreads = 1; try { realThreads = forceField.getInteger(ForceField.ForceFieldInteger.PME_REAL_THREADS); if (realThreads >= maxThreads || realThreads < 1) { throw new Exception("pme-real-threads must be < ffx.nt and greater than 0"); } concurrent = true; } catch (Exception e) { concurrent = false; } if (concurrent) { sectionThreads = 2; realSpaceThreads = realThreads; reciprocalThreads = maxThreads - realThreads; sectionTeam = new ParallelTeam(sectionThreads); realSpaceTeam = new ParallelTeam(realSpaceThreads); fftTeam = new ParallelTeam(reciprocalThreads); } else { /** * If pme-real-threads is not defined, then do real and * reciprocal space parts sequentially. */ sectionThreads = 1; realSpaceThreads = maxThreads; reciprocalThreads = maxThreads; sectionTeam = new ParallelTeam(sectionThreads); realSpaceTeam = parallelTeam; fftTeam = parallelTeam; } } realSpaceRanges = new Range[maxThreads]; initializationRegion = new InitializationRegion(maxThreads); expandInducedDipolesRegion = new ExpandInducedDipolesRegion(maxThreads); initAtomArrays(); /** * Note that we always pass on the unit cell crystal to ReciprocalSpace * instance even if the real space calculations require a * ReplicatesCrystal. */ if (aewald > 0.0 && reciprocalSpaceTerm) { reciprocalSpace = new ReciprocalSpace(this, crystal.getUnitCell(), forceField, atoms, aewald, fftTeam, parallelTeam); reciprocalEnergyRegion = new ReciprocalEnergyRegion(maxThreads); } else { reciprocalSpace = null; reciprocalEnergyRegion = null; } permanentFieldRegion = new PermanentFieldRegion(realSpaceTeam); inducedDipoleFieldRegion = new InducedDipoleFieldRegion(realSpaceTeam); directRegion = new DirectRegion(maxThreads); sorRegion = new SORRegion(maxThreads); realSpaceEnergyRegion = new RealSpaceEnergyRegion(maxThreads); reduceRegion = new ReduceRegion(maxThreads); realSpacePermTime = new long[maxThreads]; realSpaceEnergyTime = new long[maxThreads]; realSpaceSCFTime = new long[maxThreads]; /** * Generalized Kirkwood currently requires aperiodic Ewald. The GK * reaction field is added to the intra-molecular to give a * self-consistent reaction field. */ generalizedKirkwoodTerm = forceField.getBoolean(ForceFieldBoolean.GKTERM, false); if (generalizedKirkwoodTerm) { generalizedKirkwood = new GeneralizedKirkwood(forceField, atoms, this, crystal, parallelTeam); } else { generalizedKirkwood = null; } if (lambdaTerm) { StringBuilder sb = new StringBuilder(" Lambda Parameters\n"); sb.append(format(" Permanent Multipole Softcore Alpha: %5.3f\n", permLambdaAlpha)); sb.append(format(" Permanent Multipole Lambda Exponent: %5.3f\n", permLambdaExponent)); if (polarization != Polarization.NONE) { sb.append(format(" Polarization Lambda Exponent: %5.3f\n", polLambdaExponent)); sb.append( format(" Polarization Lambda Range: %5.3f .. %5.3f\n", polLambdaStart, polLambdaEnd)); sb.append(format(" Condensed SCF Without Ligand: %B\n", doNoLigandCondensedSCF)); } sb.append(format(" Vapor Electrostatics: %B\n", doLigandVaporElec)); logger.info(sb.toString()); } }
From source file:ffx.potential.nonbonded.ParticleMeshEwald.java
/** * ParticleMeshEwald constructor.//ww w . java 2 s. c om * * @param atoms the Atom array to do electrostatics on. * @param molecule the molecule number for each atom. * @param forceField the ForceField the defines the electrostatics * parameters. * @param crystal The boundary conditions. * @param elecForm The electrostatics functional form. * @param neighborList The NeighborList for both van der Waals and PME. * @param parallelTeam A ParallelTeam that delegates parallelization. */ public ParticleMeshEwald(Atom atoms[], int molecule[], ForceField forceField, Crystal crystal, NeighborList neighborList, ELEC_FORM elecForm, ParallelTeam parallelTeam) { this.atoms = atoms; this.molecule = molecule; this.forceField = forceField; this.crystal = crystal; this.parallelTeam = parallelTeam; this.neighborList = neighborList; this.elecForm = elecForm; neighborLists = neighborList.getNeighborList(); permanentSchedule = neighborList.getPairwiseSchedule(); nAtoms = atoms.length; nSymm = crystal.spaceGroup.getNumberOfSymOps(); maxThreads = parallelTeam.getThreadCount(); polsor = forceField.getDouble(ForceFieldDouble.POLAR_SOR, 0.70); poleps = forceField.getDouble(ForceFieldDouble.POLAR_EPS, 1e-5); if (elecForm == ELEC_FORM.PAM) { m12scale = forceField.getDouble(ForceFieldDouble.MPOLE_12_SCALE, 0.0); m13scale = forceField.getDouble(ForceFieldDouble.MPOLE_13_SCALE, 0.0); m14scale = forceField.getDouble(ForceFieldDouble.MPOLE_14_SCALE, 0.4); m15scale = forceField.getDouble(ForceFieldDouble.MPOLE_15_SCALE, 0.8); } else { m12scale = forceField.getDouble(ForceFieldDouble.MPOLE_12_SCALE, 0.0); m13scale = forceField.getDouble(ForceFieldDouble.MPOLE_13_SCALE, 0.0); m14scale = forceField.getDouble(ForceFieldDouble.MPOLE_14_SCALE, 0.5); m15scale = forceField.getDouble(ForceFieldDouble.MPOLE_15_SCALE, 1.0); } d11scale = forceField.getDouble(ForceFieldDouble.DIRECT_11_SCALE, 0.0); p12scale = forceField.getDouble(ForceFieldDouble.POLAR_12_SCALE, 0.0); p13scale = forceField.getDouble(ForceFieldDouble.POLAR_13_SCALE, 0.0); useCharges = forceField.getBoolean(ForceFieldBoolean.USE_CHARGES, true); useDipoles = forceField.getBoolean(ForceFieldBoolean.USE_DIPOLES, true); useQuadrupoles = forceField.getBoolean(ForceFieldBoolean.USE_QUADRUPOLES, true); rotateMultipoles = forceField.getBoolean(ForceFieldBoolean.ROTATE_MULTIPOLES, true); lambdaTerm = forceField.getBoolean(ForceFieldBoolean.LAMBDATERM, false); if (!crystal.aperiodic()) { off = forceField.getDouble(ForceFieldDouble.EWALD_CUTOFF, 7.0); } else { off = forceField.getDouble(ForceFieldDouble.EWALD_CUTOFF, 100.0); } double ewald_precision = forceField.getDouble(ForceFieldDouble.EWALD_PRECISION, 1.0e-8); aewald = forceField.getDouble(ForceFieldDouble.EWALD_ALPHA, ewaldCoefficient(off, ewald_precision)); setEwaldParameters(off, aewald); String predictor = forceField.getString(ForceFieldString.SCF_PREDICTOR, "NONE"); try { predictor = predictor.replaceAll("-", "_").toUpperCase(); scfPredictor = SCFPredictor.valueOf(predictor); } catch (Exception e) { scfPredictor = SCFPredictor.NONE; } if (scfPredictor != SCFPredictor.NONE) { predictorCount = 0; int defaultOrder = 6; predictorOrder = forceField.getInteger(ForceFieldInteger.SCF_PREDICTOR_ORDER, defaultOrder); if (scfPredictor == SCFPredictor.LS) { leastSquaresPredictor = new LeastSquaresPredictor(); double eps = 1.0e-4; leastSquaresOptimizer = new LevenbergMarquardtOptimizer(new SimpleVectorValueChecker(eps, eps)); } else if (scfPredictor == SCFPredictor.ASPC) { predictorOrder = 6; } predictorStartIndex = 0; } String algorithm = forceField.getString(ForceFieldString.SCF_ALGORITHM, "CG"); try { algorithm = algorithm.replaceAll("-", "_").toUpperCase(); scfAlgorithm = SCFAlgorithm.valueOf(algorithm); } catch (Exception e) { scfAlgorithm = SCFAlgorithm.CG; } /** * The size of the preconditioner neighbor list depends on the size of * the preconditioner cutoff. */ if (scfAlgorithm == SCFAlgorithm.CG) { inducedDipolePreconditionerRegion = new InducedDipolePreconditionerRegion(maxThreads); pcgRegion = new PCGRegion(maxThreads); pcgInitRegion1 = new PCGInitRegion1(maxThreads); pcgInitRegion2 = new PCGInitRegion2(maxThreads); pcgIterRegion1 = new PCGIterRegion1(maxThreads); pcgIterRegion2 = new PCGIterRegion2(maxThreads); boolean preconditioner = forceField.getBoolean(ForceFieldBoolean.USE_SCF_PRECONDITIONER, true); if (preconditioner) { preconditionerCutoff = forceField.getDouble(ForceFieldDouble.CG_PRECONDITIONER_CUTOFF, 4.5); preconditionerEwald = forceField.getDouble(ForceFieldDouble.CG_PRECONDITIONER_EWALD, 0.0); } else { preconditionerCutoff = 0.0; } } else { preconditionerCutoff = 0.0; inducedDipolePreconditionerRegion = null; pcgRegion = null; pcgInitRegion1 = null; pcgInitRegion2 = null; pcgIterRegion1 = null; pcgIterRegion2 = null; } if (lambdaTerm) { /** * Values of PERMANENT_LAMBDA_ALPHA below 2 can lead to unstable * trajectories. */ permLambdaAlpha = forceField.getDouble(ForceFieldDouble.PERMANENT_LAMBDA_ALPHA, 2.0); if (permLambdaAlpha < 0.0 || permLambdaAlpha > 3.0) { permLambdaAlpha = 2.0; } /** * A PERMANENT_LAMBDA_EXPONENT of 1 gives linear charging of the * permanent electrostatics, which is most efficient. A quadratic * schedule (PERMANENT_LAMBDA_EXPONENT) also works, but the dU/dL * forces near lambda=1 are may be larger by a factor of 2. */ permLambdaExponent = forceField.getDouble(ForceFieldDouble.PERMANENT_LAMBDA_EXPONENT, 1.0); if (permLambdaExponent < 1.0) { permLambdaExponent = 1.0; } /** * A POLARIZATION_LAMBDA_EXPONENT of 2 gives a non-zero d2U/dL2 at * the beginning of the polarization schedule. Choosing a power of 3 * or greater ensures a smooth dU/dL and d2U/dL2 over the schedule. */ polLambdaExponent = forceField.getDouble(ForceFieldDouble.POLARIZATION_LAMBDA_EXPONENT, 3.0); if (polLambdaExponent < 0.0) { polLambdaExponent = 0.0; } /** * The POLARIZATION_LAMBDA_START defines the point in the lambda * schedule when the condensed phase polarization of the ligand * begins to be turned on. If the condensed phase polarization is * considered near lambda=0, then SCF convergence is slow, even with * Thole damping. In addition, 2 (instead of 1) condensed phase SCF * calculations are necessary from the beginning of the window to * lambda=1. */ polLambdaStart = forceField.getDouble(ForceFieldDouble.POLARIZATION_LAMBDA_START, 0.75); if (polLambdaStart < 0.0 || polLambdaStart > 0.9) { polLambdaStart = 0.75; } /** * The POLARIZATION_LAMBDA_END defines the point in the lambda * schedule when the condensed phase polarization of ligand has been * completely turned on. Values other than 1.0 have not been tested. */ polLambdaEnd = forceField.getDouble(ForceFieldDouble.POLARIZATION_LAMBDA_END, 1.0); if (polLambdaEnd < polLambdaStart || polLambdaEnd > 1.0 || polLambdaEnd - polLambdaStart < 0.3) { polLambdaEnd = 1.0; } /** * The LAMBDA_VAPOR_ELEC defines if intramolecular electrostatics of * the ligand in vapor will be considered. */ doLigandVaporElec = forceField.getBoolean(ForceFieldBoolean.LIGAND_VAPOR_ELEC, true); doNoLigandCondensedSCF = forceField.getBoolean(ForceFieldBoolean.NO_LIGAND_CONDENSED_SCF, true); /** * Flag to indicate application of an intermolecular softcore * potential. */ intermolecularSoftcore = forceField.getBoolean(ForceField.ForceFieldBoolean.INTERMOLECULAR_SOFTCORE, false); } String polar = forceField.getString(ForceFieldString.POLARIZATION, "MUTUAL"); if (elecForm == ELEC_FORM.FIXED_CHARGE) { polar = "NONE"; } boolean polarizationTerm = forceField.getBoolean(ForceFieldBoolean.POLARIZETERM, true); if (polarizationTerm == false || polar.equalsIgnoreCase("NONE")) { polarization = Polarization.NONE; } else if (polar.equalsIgnoreCase("DIRECT")) { polarization = Polarization.DIRECT; } else { polarization = Polarization.MUTUAL; } cudaFFT = forceField.getBoolean(ForceField.ForceFieldBoolean.CUDAFFT, false); if (lambdaTerm) { shareddEdLambda = new SharedDouble(); sharedd2EdLambda2 = new SharedDouble(); } else { shareddEdLambda = null; sharedd2EdLambda2 = null; lambdaGrad = null; lambdaTorque = null; vaporCrystal = null; vaporLists = null; vaporPermanentSchedule = null; vaporEwaldSchedule = null; vacuumRanges = null; } if (logger.isLoggable(Level.INFO)) { StringBuilder sb = new StringBuilder("\n Electrostatics\n"); sb.append(format(" Polarization: %8s\n", polarization.toString())); if (polarization == Polarization.MUTUAL) { sb.append(format(" SCF Convergence Criteria: %8.3e\n", poleps)); sb.append(format(" SCF Predictor: %8s\n", scfPredictor)); sb.append(format(" SCF Algorithm: %8s\n", scfAlgorithm)); if (scfAlgorithm == SCFAlgorithm.SOR) { sb.append(format(" SOR Parameter: %8.3f\n", polsor)); } else { sb.append(format(" CG Preconditioner Cut-Off: %8.3f\n", preconditionerCutoff)); sb.append(format(" CG Preconditioner Ewald Coefficient:%8.3f\n", preconditionerEwald)); } } if (aewald > 0.0) { sb.append(" Particle-mesh Ewald\n"); sb.append(format(" Ewald Coefficient: %8.3f\n", aewald)); sb.append(format(" Particle Cut-Off: %8.3f (A)", off)); } else { sb.append(format(" Electrostatics Cut-Off: %8.3f (A)\n", off)); } logger.info(sb.toString()); } if (cudaFFT) { sectionThreads = 2; realSpaceThreads = parallelTeam.getThreadCount(); reciprocalThreads = 1; sectionTeam = new ParallelTeam(sectionThreads); realSpaceTeam = parallelTeam; fftTeam = new ParallelTeam(reciprocalThreads); } else { boolean concurrent; int realThreads = 1; try { realThreads = forceField.getInteger(ForceField.ForceFieldInteger.PME_REAL_THREADS); if (realThreads >= maxThreads || realThreads < 1) { throw new Exception("pme-real-threads must be < ffx.nt and greater than 0"); } concurrent = true; } catch (Exception e) { concurrent = false; } if (concurrent) { sectionThreads = 2; realSpaceThreads = realThreads; reciprocalThreads = maxThreads - realThreads; sectionTeam = new ParallelTeam(sectionThreads); realSpaceTeam = new ParallelTeam(realSpaceThreads); fftTeam = new ParallelTeam(reciprocalThreads); } else { /** * If pme-real-threads is not defined, then do real and * reciprocal space parts sequentially. */ sectionThreads = 1; realSpaceThreads = maxThreads; reciprocalThreads = maxThreads; sectionTeam = new ParallelTeam(sectionThreads); realSpaceTeam = parallelTeam; fftTeam = parallelTeam; } } realSpaceRanges = new Range[maxThreads]; initializationRegion = new InitializationRegion(maxThreads); expandInducedDipolesRegion = new ExpandInducedDipolesRegion(maxThreads); initAtomArrays(); /** * Note that we always pass on the unit cell crystal to ReciprocalSpace * instance even if the real space calculations require a * ReplicatesCrystal. */ if (aewald > 0.0) { reciprocalSpace = new ReciprocalSpace(this, crystal.getUnitCell(), forceField, atoms, aewald, fftTeam, parallelTeam); reciprocalEnergyRegion = new ReciprocalEnergyRegion(maxThreads); } else { reciprocalSpace = null; reciprocalEnergyRegion = null; } permanentFieldRegion = new PermanentFieldRegion(realSpaceTeam); inducedDipoleFieldRegion = new InducedDipoleFieldRegion(realSpaceTeam); directRegion = new DirectRegion(maxThreads); sorRegion = new SORRegion(maxThreads); realSpaceEnergyRegion = new RealSpaceEnergyRegion(maxThreads); reduceRegion = new ReduceRegion(maxThreads); realSpacePermTime = new long[maxThreads]; realSpaceEnergyTime = new long[maxThreads]; realSpaceSCFTime = new long[maxThreads]; /** * Generalized Kirkwood currently requires aperiodic Ewald. The GK * reaction field is added to the intra-molecular to give a * self-consistent reaction field. */ generalizedKirkwoodTerm = forceField.getBoolean(ForceFieldBoolean.GKTERM, false); if (generalizedKirkwoodTerm) { generalizedKirkwood = new GeneralizedKirkwood(forceField, atoms, this, crystal, parallelTeam); } else { generalizedKirkwood = null; } if (lambdaTerm) { StringBuilder sb = new StringBuilder(" Lambda Parameters\n"); sb.append(format(" Permanent Multipole Softcore Alpha: %5.3f\n", permLambdaAlpha)); sb.append(format(" Permanent Multipole Lambda Exponent: %5.3f\n", permLambdaExponent)); if (polarization != Polarization.NONE) { sb.append(format(" Polarization Lambda Exponent: %5.3f\n", polLambdaExponent)); sb.append( format(" Polarization Lambda Range: %5.3f .. %5.3f\n", polLambdaStart, polLambdaEnd)); sb.append(format(" Condensed SCF Without Ligand: %B\n", doNoLigandCondensedSCF)); } sb.append(format(" Vapor Electrostatics: %B\n", doLigandVaporElec)); logger.info(sb.toString()); } }
From source file:ffx.potential.nonbonded.ParticleMeshEwaldQI.java
/** * ParticleMeshEwald constructor.//from w ww . j av a 2s .co m * * @param atoms the Atom array to do electrostatics on. * @param molecule the molecule number for each atom. * @param forceField the ForceField the defines the electrostatics * parameters. * @param crystal The boundary conditions. * @param elecForm The electrostatics functional form. * @param neighborList The NeighborList for both van der Waals and PME. * @param parallelTeam A ParallelTeam that delegates parallelization. */ public ParticleMeshEwaldQI(Atom atoms[], int molecule[], ForceField forceField, Crystal crystal, NeighborList neighborList, ELEC_FORM elecForm, ParallelTeam parallelTeam) { /* REM Used to require the dlAlphaMode == FACTORED. ie. dlAlpha /= -2.0, d2lAlpha /= -2.0; */ bufferCoords = COORDINATES.QI; this.atoms = atoms; this.molecule = molecule; this.forceField = forceField; this.crystal = crystal; this.parallelTeam = parallelTeam; this.neighborList = neighborList; this.elecForm = elecForm; neighborLists = neighborList.getNeighborList(); permanentSchedule = neighborList.getPairwiseSchedule(); nAtoms = atoms.length; nSymm = crystal.spaceGroup.getNumberOfSymOps(); maxThreads = parallelTeam.getThreadCount() + 1; polsor = forceField.getDouble(ForceFieldDouble.POLAR_SOR, 0.70); poleps = forceField.getDouble(ForceFieldDouble.POLAR_EPS, 1e-5); if (elecForm == ELEC_FORM.PAM) { m12scale = forceField.getDouble(ForceFieldDouble.MPOLE_12_SCALE, 0.0); m13scale = forceField.getDouble(ForceFieldDouble.MPOLE_13_SCALE, 0.0); m14scale = forceField.getDouble(ForceFieldDouble.MPOLE_14_SCALE, 0.4); m15scale = forceField.getDouble(ForceFieldDouble.MPOLE_15_SCALE, 0.8); } else { double mpole14 = 0.5; String name = forceField.toString().toUpperCase(); if (name.contains("AMBER")) { mpole14 = 1.0 / 1.2; } m12scale = forceField.getDouble(ForceFieldDouble.MPOLE_12_SCALE, 0.0); m13scale = forceField.getDouble(ForceFieldDouble.MPOLE_13_SCALE, 0.0); m14scale = forceField.getDouble(ForceFieldDouble.MPOLE_14_SCALE, mpole14); m15scale = forceField.getDouble(ForceFieldDouble.MPOLE_15_SCALE, 1.0); } d11scale = forceField.getDouble(ForceFieldDouble.DIRECT_11_SCALE, 0.0); p12scale = forceField.getDouble(ForceFieldDouble.POLAR_12_SCALE, 0.0); p13scale = forceField.getDouble(ForceFieldDouble.POLAR_13_SCALE, 0.0); useCharges = forceField.getBoolean(ForceFieldBoolean.USE_CHARGES, true); useDipoles = forceField.getBoolean(ForceFieldBoolean.USE_DIPOLES, true); useQuadrupoles = forceField.getBoolean(ForceFieldBoolean.USE_QUADRUPOLES, true); rotateMultipoles = forceField.getBoolean(ForceFieldBoolean.ROTATE_MULTIPOLES, true); lambdaTerm = forceField.getBoolean(ForceFieldBoolean.LAMBDATERM, false); if (!crystal.aperiodic()) { off = forceField.getDouble(ForceFieldDouble.EWALD_CUTOFF, 7.0); } else { off = forceField.getDouble(ForceFieldDouble.EWALD_CUTOFF, 1000.0); } double ewaldPrecision = forceField.getDouble(ForceFieldDouble.EWALD_PRECISION, 1.0e-8); aewald = forceField.getDouble(ForceFieldDouble.EWALD_ALPHA, ewaldCoefficient(off, ewaldPrecision)); setEwaldParameters(off, aewald); reciprocalSpaceTerm = forceField.getBoolean(ForceFieldBoolean.RECIPTERM, true); String predictor = forceField.getString(ForceFieldString.SCF_PREDICTOR, "NONE"); try { predictor = predictor.replaceAll("-", "_").toUpperCase(); scfPredictor = SCFPredictor.valueOf(predictor); } catch (Exception e) { scfPredictor = SCFPredictor.NONE; } if (scfPredictor != SCFPredictor.NONE) { predictorCount = 0; int defaultOrder = 6; predictorOrder = forceField.getInteger(ForceFieldInteger.SCF_PREDICTOR_ORDER, defaultOrder); if (scfPredictor == SCFPredictor.LS) { leastSquaresPredictor = new LeastSquaresPredictor(); double eps = 1.0e-4; leastSquaresOptimizer = new LevenbergMarquardtOptimizer(new SimpleVectorValueChecker(eps, eps)); } else if (scfPredictor == SCFPredictor.ASPC) { predictorOrder = 6; } predictorStartIndex = 0; } String algorithm = forceField.getString(ForceFieldString.SCF_ALGORITHM, "CG"); try { algorithm = algorithm.replaceAll("-", "_").toUpperCase(); scfAlgorithm = SCFAlgorithm.valueOf(algorithm); } catch (Exception e) { scfAlgorithm = SCFAlgorithm.CG; } /** * The size of the preconditioner neighbor list depends on the size of * the preconditioner cutoff. */ if (scfAlgorithm == SCFAlgorithm.CG) { inducedDipolePreconditionerRegion = new InducedDipolePreconditionerRegion(maxThreads); pcgRegion = new PCGRegion(maxThreads); pcgInitRegion1 = new PCGInitRegion1(maxThreads); pcgInitRegion2 = new PCGInitRegion2(maxThreads); pcgIterRegion1 = new PCGIterRegion1(maxThreads); pcgIterRegion2 = new PCGIterRegion2(maxThreads); boolean preconditioner = forceField.getBoolean(ForceFieldBoolean.USE_SCF_PRECONDITIONER, true); if (preconditioner) { preconditionerCutoff = forceField.getDouble(ForceFieldDouble.CG_PRECONDITIONER_CUTOFF, 4.5); preconditionerEwald = forceField.getDouble(ForceFieldDouble.CG_PRECONDITIONER_EWALD, 0.0); } else { preconditionerCutoff = 0.0; } } else { preconditionerCutoff = 0.0; inducedDipolePreconditionerRegion = null; pcgRegion = null; pcgInitRegion1 = null; pcgInitRegion2 = null; pcgIterRegion1 = null; pcgIterRegion2 = null; } if (lambdaTerm || esvTerm) { /** * Values of PERMANENT_LAMBDA_ALPHA below 2 can lead to unstable * trajectories. */ permLambdaAlpha = forceField.getDouble(ForceFieldDouble.PERMANENT_LAMBDA_ALPHA, 2.0); if (permLambdaAlpha < 0.0 || permLambdaAlpha > 3.0) { permLambdaAlpha = 2.0; } /** * A PERMANENT_LAMBDA_EXPONENT of 1 gives linear charging of the * permanent electrostatics, which is most efficient. A quadratic * schedule (PERMANENT_LAMBDA_EXPONENT) also works, but the dU/dL * forces near lambda=1 are may be larger by a factor of 2. */ permLambdaExponent = forceField.getDouble(ForceFieldDouble.PERMANENT_LAMBDA_EXPONENT, 1.0); if (permLambdaExponent < 1.0) { permLambdaExponent = 1.0; } /** * A POLARIZATION_LAMBDA_EXPONENT of 2 gives a non-zero d2U/dL2 at * the beginning of the polarization schedule. Choosing a power of 3 * or greater ensures a smooth dU/dL and d2U/dL2 over the schedule. */ polLambdaExponent = forceField.getDouble(ForceFieldDouble.POLARIZATION_LAMBDA_EXPONENT, 3.0); if (polLambdaExponent < 0.0) { polLambdaExponent = 0.0; } /** * The POLARIZATION_LAMBDA_START defines the point in the lambda * schedule when the condensed phase polarization of the ligand * begins to be turned on. If the condensed phase polarization is * considered near lambda=0, then SCF convergence is slow, even with * Thole damping. In addition, 2 (instead of 1) condensed phase SCF * calculations are necessary from the beginning of the window to * lambda=1. */ polLambdaStart = forceField.getDouble(ForceFieldDouble.POLARIZATION_LAMBDA_START, 0.75); if (polLambdaStart < 0.0 || polLambdaStart > 0.9) { polLambdaStart = 0.75; } /** * The POLARIZATION_LAMBDA_END defines the point in the lambda * schedule when the condensed phase polarization of ligand has been * completely turned on. Values other than 1.0 have not been tested. */ polLambdaEnd = forceField.getDouble(ForceFieldDouble.POLARIZATION_LAMBDA_END, 1.0); if (polLambdaEnd < polLambdaStart || polLambdaEnd > 1.0 || polLambdaEnd - polLambdaStart < 0.3) { polLambdaEnd = 1.0; } /** * The LAMBDA_VAPOR_ELEC defines if intramolecular electrostatics of * the ligand in vapor will be considered. */ doLigandVaporElec = forceField.getBoolean(ForceFieldBoolean.LIGAND_VAPOR_ELEC, true); doNoLigandCondensedSCF = forceField.getBoolean(ForceFieldBoolean.NO_LIGAND_CONDENSED_SCF, true); /** * Flag to indicate application of an intermolecular softcore * potential. */ intermolecularSoftcore = forceField.getBoolean(ForceField.ForceFieldBoolean.INTERMOLECULAR_SOFTCORE, false); intramolecularSoftcore = forceField.getBoolean(ForceField.ForceFieldBoolean.INTRAMOLECULAR_SOFTCORE, false); } String polar = forceField.getString(ForceFieldString.POLARIZATION, "MUTUAL"); if (elecForm == ELEC_FORM.FIXED_CHARGE) { polar = "NONE"; } boolean polarizationTerm = forceField.getBoolean(ForceFieldBoolean.POLARIZETERM, true); if (polarizationTerm == false || polar.equalsIgnoreCase("NONE")) { polarization = Polarization.NONE; } else if (polar.equalsIgnoreCase("DIRECT")) { polarization = Polarization.DIRECT; } else { polarization = Polarization.MUTUAL; } String temp = forceField.getString(ForceField.ForceFieldString.FFT_METHOD, "PJ"); FFTMethod method; try { method = ReciprocalSpace.FFTMethod.valueOf(temp.toUpperCase().trim()); } catch (Exception e) { method = ReciprocalSpace.FFTMethod.PJ; } gpuFFT = method != FFTMethod.PJ; if (lambdaTerm) { shareddEdLambda = new SharedDouble(); sharedd2EdLambda2 = new SharedDouble(); } else { shareddEdLambda = null; sharedd2EdLambda2 = null; lambdaGrad = null; lambdaTorque = null; vaporCrystal = null; vaporLists = null; vaporPermanentSchedule = null; vaporEwaldSchedule = null; vacuumRanges = null; } if (logger.isLoggable(Level.INFO)) { StringBuilder sb = new StringBuilder("\n Electrostatics\n"); sb.append(format(" Polarization: %8s\n", polarization.toString())); if (polarization == Polarization.MUTUAL) { sb.append(format(" SCF Convergence Criteria: %8.3e\n", poleps)); sb.append(format(" SCF Predictor: %8s\n", scfPredictor)); sb.append(format(" SCF Algorithm: %8s\n", scfAlgorithm)); if (scfAlgorithm == SCFAlgorithm.SOR) { sb.append(format(" SOR Parameter: %8.3f\n", polsor)); } else { sb.append(format(" CG Preconditioner Cut-Off: %8.3f\n", preconditionerCutoff)); sb.append(format(" CG Preconditioner Ewald Coefficient:%8.3f\n", preconditionerEwald)); } } if (aewald > 0.0) { sb.append(" Particle-mesh Ewald\n"); sb.append(format(" Ewald Coefficient: %8.3f\n", aewald)); sb.append(format(" Particle Cut-Off: %8.3f (A)", off)); } else { sb.append(format(" Electrostatics Cut-Off: %8.3f (A)\n", off)); } logger.info(sb.toString()); } StringBuilder config = new StringBuilder(); config.append(format("\n Quasi-Internal PME Settings\n")); config.append(format(" Debug,Verbose,dbgI&K: %5b %5b %5d %5d\n", DEBUG(), VERBOSE(), debugIntI().orElse(-1), debugIntK().orElse(-1))); config.append(format(" Buffer Coords: %5s\n", bufferCoords.toString())); config.append(format(" Chrg,Dipl,Quad: %5b %5b %5b\n", useCharges, useDipoles, useQuadrupoles)); logger.info(config.toString()); if (gpuFFT) { sectionThreads = 2; realSpaceThreads = parallelTeam.getThreadCount(); reciprocalThreads = 1; sectionTeam = new ParallelTeam(sectionThreads); realSpaceTeam = parallelTeam; fftTeam = new ParallelTeam(reciprocalThreads); } else { boolean concurrent; int realThreads = 1; try { realThreads = forceField.getInteger(ForceField.ForceFieldInteger.PME_REAL_THREADS); if (realThreads >= maxThreads || realThreads < 1) { throw new Exception("pme-real-threads must be < ffx.nt and greater than 0"); } concurrent = true; } catch (Exception e) { concurrent = false; } if (concurrent) { sectionThreads = 2; realSpaceThreads = realThreads; reciprocalThreads = maxThreads - realThreads; sectionTeam = new ParallelTeam(sectionThreads); realSpaceTeam = new ParallelTeam(realSpaceThreads); fftTeam = new ParallelTeam(reciprocalThreads); } else { /** * If pme-real-threads is not defined, then do real and * reciprocal space parts sequentially. */ sectionThreads = 1; realSpaceThreads = maxThreads; reciprocalThreads = maxThreads; sectionTeam = new ParallelTeam(sectionThreads); realSpaceTeam = parallelTeam; fftTeam = parallelTeam; } } realSpaceRanges = new Range[maxThreads]; initializationRegion = new InitializationRegion(maxThreads); expandInducedDipolesRegion = new ExpandInducedDipolesRegion(maxThreads); initAtomArrays(); /** * Note that we always pass on the unit cell crystal to ReciprocalSpace * instance even if the real space calculations require a * ReplicatesCrystal. */ if (aewald > 0.0 && reciprocalSpaceTerm) { reciprocalSpace = new ReciprocalSpace(this, crystal.getUnitCell(), forceField, atoms, aewald, fftTeam, parallelTeam); reciprocalEnergyRegion = new ReciprocalEnergyRegion(maxThreads); } else { reciprocalSpace = null; reciprocalEnergyRegion = null; } permanentFieldRegion = new PermanentFieldRegion(realSpaceTeam); inducedDipoleFieldRegion = new InducedDipoleFieldRegion(realSpaceTeam); directRegion = new DirectRegion(maxThreads); sorRegion = new SORRegion(maxThreads); realSpaceEnergyRegionQI = new RealSpaceEnergyRegionQI(maxThreads); reduceRegion = new ReduceRegion(maxThreads); realSpacePermTime = new long[maxThreads]; realSpaceEnergyTime = new long[maxThreads]; realSpaceSCFTime = new long[maxThreads]; /** * Generalized Kirkwood currently requires aperiodic Ewald. The GK * reaction field is added to the intra-molecular to give a * self-consistent reaction field. */ generalizedKirkwoodTerm = forceField.getBoolean(ForceFieldBoolean.GKTERM, false); if (generalizedKirkwoodTerm) { generalizedKirkwood = new GeneralizedKirkwood(forceField, atoms, this, crystal, parallelTeam); } else { generalizedKirkwood = null; } if (lambdaTerm) { StringBuilder sb = new StringBuilder(" Lambda Parameters\n"); sb.append(format(" Permanent Multipole Softcore Alpha: %5.3f\n", permLambdaAlpha)); sb.append(format(" Permanent Multipole Lambda Exponent: %5.3f\n", permLambdaExponent)); if (polarization != Polarization.NONE) { sb.append(format(" Polarization Lambda Exponent: %5.3f\n", polLambdaExponent)); sb.append( format(" Polarization Lambda Range: %5.3f .. %5.3f\n", polLambdaStart, polLambdaEnd)); sb.append(format(" Condensed SCF Without Ligand: %B\n", doNoLigandCondensedSCF)); } sb.append(format(" Vapor Electrostatics: %B\n", doLigandVaporElec)); logger.info(sb.toString()); } if (esvTerm) { StringBuilder sb = new StringBuilder(" ESV Parameters\n"); sb.append(format(" Permanent Multipole Softcore Alpha: %5.3f\n", permLambdaAlpha)); sb.append(format(" Permanent Multipole lambda exponent constrained to unity.\n")); if (polarization != Polarization.NONE) { throw new UnsupportedOperationException(); } logger.info(sb.toString()); } }