From feb7daea21ec62ad8f3e682dce60a9d6f6fdfb9b Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Thu, 23 Jan 2014 14:51:48 -0800 Subject: [PATCH 01/34] Fix to polarization correction --- src/forcefield.py | 3 ++- src/liquid.py | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/forcefield.py b/src/forcefield.py index 9953b7cf9..8691ee31c 100644 --- a/src/forcefield.py +++ b/src/forcefield.py @@ -1113,7 +1113,8 @@ def build_qtrans2(tq, qid, qmap): def list_map(self): """ Create the plist, which is like a reversed version of the parameter map. More convenient for printing. """ if len(self.map) == 0: - warn_press_key('The parameter map has no elements (Okay if we are not actually tuning any parameters.)') + #warn_press_key('The parameter map has no elements (Okay if we are not actually tuning any parameters.)') + logger.warning('The parameter map has no elements (Okay if we are not actually tuning any parameters.)\n') else: self.plist = [[] for j in range(max([self.map[i] for i in self.map])+1)] for i in self.map: diff --git a/src/liquid.py b/src/liquid.py index 2495c0178..1f5a920df 100644 --- a/src/liquid.py +++ b/src/liquid.py @@ -257,6 +257,7 @@ def npt_simulation(self, temperature, pressure, simnum): output_files = ['npt_result.p.bz2', 'npt.out'] + self.extra_output, tgt=self) def polarization_correction(self,mvals): + self.FF.make(mvals) d = self.gas_engine.multipole_moments(optimize=True)['dipole'] if not in_fd(): logger.info("The molecular dipole moment is % .3f debye\n" % np.linalg.norm(d)) From bc3a415a78bd46db55e1a4904170f62e36c20a30 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Tue, 11 Feb 2014 14:54:18 -0800 Subject: [PATCH 02/34] Begin adding interface for RP contraction --- src/data/npt.py | 6 +++--- src/forcefield.py | 2 +- src/openmmio.py | 52 +++++++++++++++++++++++++++++++++++++++-------- src/parser.py | 4 ++-- 4 files changed, 49 insertions(+), 15 deletions(-) diff --git a/src/data/npt.py b/src/data/npt.py index 659f44c69..a747a671f 100755 --- a/src/data/npt.py +++ b/src/data/npt.py @@ -311,7 +311,7 @@ def main(): # Number of threads, multiple timestep integrator, anisotropic box etc. threads = TgtOptions.get('md_threads', 1) mts = TgtOptions.get('mts_integrator', 0) - rpmd_beads = TgtOptions.get('rpmd_beads', 0) + rpmd_beads = TgtOptions.get('rpmd_beads', []) force_cuda = TgtOptions.get('force_cuda', 0) anisotropic = TgtOptions.get('anisotropic_box', 0) minimize = TgtOptions.get('minimize_energy', 1) @@ -365,7 +365,7 @@ def main(): EngOpts["gas"]["gmx_top"] = os.path.splitext(gas_fnm)[0] + ".top" EngOpts["gas"]["gmx_mdp"] = os.path.splitext(gas_fnm)[0] + ".mdp" if force_cuda: logger.warn("force_cuda option has no effect on Gromacs engine.") - if rpmd_beads > 0: raise RuntimeError("Gromacs cannot handle RPMD.") + if len(rpmd_beads) > 0: raise RuntimeError("Gromacs cannot handle RPMD.") if mts: logger.warn("Gromacs not configured for multiple timestep integrator.") if anisotropic: logger.warn("Gromacs not configured for anisotropic box scaling.") elif engname == "tinker": @@ -374,7 +374,7 @@ def main(): EngOpts["liquid"]["tinker_key"] = os.path.splitext(liquid_fnm)[0] + ".key" EngOpts["gas"]["tinker_key"] = os.path.splitext(gas_fnm)[0] + ".key" if force_cuda: logger.warn("force_cuda option has no effect on Tinker engine.") - if rpmd_beads > 0: raise RuntimeError("TINKER cannot handle RPMD.") + if len(rpmd_beads) > 0: raise RuntimeError("TINKER cannot handle RPMD.") if mts: logger.warn("Tinker not configured for multiple timestep integrator.") EngOpts["liquid"].update(GenOpts) EngOpts["gas"].update(GenOpts) diff --git a/src/forcefield.py b/src/forcefield.py index af42ab598..4ab7ae62a 100644 --- a/src/forcefield.py +++ b/src/forcefield.py @@ -223,7 +223,7 @@ def __init__(self, options, verbose=True): ## AMOEBA mutual dipole convergence tolerance. self.set_option(options, 'amoeba_eps') ## Switch for rigid water molecules - self.set_option(options, 'rigid_water') + self.set_option(options, 'rigid_water', forceprint=True) ## Bypass the transformation and use physical parameters directly self.set_option(options, 'use_pvals') diff --git a/src/openmmio.py b/src/openmmio.py index cd95486be..0d4129d0e 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -566,7 +566,7 @@ def prepare(self, pbc=False, mmopts={}, **kwargs): self.AtomLists['ResidueNumber'] = [a.residue.index for a in Atoms] self.AtomMask = [a == 'A' for a in self.AtomLists['ParticleType']] - def create_simulation(self, timestep=1.0, faststep=0.25, temperature=None, pressure=None, anisotropic=False, mts=False, collision=1.0, nbarostat=25, rpmd_beads=0, **kwargs): + def create_simulation(self, timestep=1.0, faststep=0.25, temperature=None, pressure=None, anisotropic=False, mts=False, collision=1.0, nbarostat=25, rpmd_beads=[], **kwargs): """ Create simulation object. Note that this also takes in some @@ -581,20 +581,54 @@ def create_simulation(self, timestep=1.0, faststep=0.25, temperature=None, press if temperature: ## If temperature control is turned on, then run Langevin dynamics. if mts: - if rpmd_beads > 0: + if len(rpmd_beads) > 0: raise RuntimeError("No multiple timestep integrator without temperature control.") integrator = MTSVVVRIntegrator(temperature*kelvin, collision/picosecond, timestep*femtosecond, self.system, ninnersteps=int(timestep/faststep)) else: - if rpmd_beads > 0: - logger.info("Creating RPMD integrator with %i beads.\n" % rpmd_beads) - self.tdiv = rpmd_beads - integrator = RPMDIntegrator(rpmd_beads, temperature*kelvin, collision/picosecond, timestep*femtosecond) + if len(rpmd_beads) > 0: + self.tdiv = rpmd_beads[0] + if len(rpmd_beads) == 1: + logger.info("Creating RPMD integrator with %i beads.\n" % rpmd_beads[0]) + integrator = RPMDIntegrator(rpmd_beads[0], temperature*kelvin, collision/picosecond, timestep*femtosecond) + elif len(rpmd_beads) == 2: + contract = False + for frc in self.system.getForces(): + if any([isinstance(frc, fc) for fc in [NonbondedForce, AmoebaMultipoleForce, AmoebaVdwForce, CustomNonbondedForce]]): + contract = True + frc.setForceGroup(1) + if contract: + logger.info("Creating RPMD integrator with %i beads (NB forces contracted to %i).\n" % (rpmd_beads[0], rpmd_beads[1])) + integrator = RPMDIntegrator(rpmd_beads[0], temperature*kelvin, collision/picosecond, timestep*femtosecond, {1:rpmd_beads[1]}) + else: + logger.info("Creating RPMD integrator with %i beads (no NB forces to contract).\n" % (rpmd_beads[0])) + integrator = RPMDIntegrator(rpmd_beads[0], temperature*kelvin, collision/picosecond, timestep*femtosecond) + elif len(rpmd_beads) == 3: + contract = False + contract_recip = False + for frc in self.system.getForces(): + if any([isinstance(frc, fc) for fc in [NonbondedForce, AmoebaMultipoleForce, AmoebaVdwForce, CustomNonbondedForce]]): + contract = True + frc.setForceGroup(1) + if isinstance(frc, NonbondedForce): + contract_recip = True + frc.setReciprocalSpaceForceGroup(2) + if contract_recip: + logger.info("Creating RPMD integrator with %i beads (NB/Recip forces contracted to %i/%i).\n" % (rpmd_beads[0], rpmd_beads[1], rpmd_beads[2])) + integrator = RPMDIntegrator(rpmd_beads[0], temperature*kelvin, collision/picosecond, timestep*femtosecond, {1:rpmd_beads[1], 2:rpmd_beads[2]}) + elif contract: + logger.info("Creating RPMD integrator with %i beads (NB forces contracted to %i, no Recip).\n" % (rpmd_beads[0], rpmd_beads[1])) + integrator = RPMDIntegrator(rpmd_beads[0], temperature*kelvin, collision/picosecond, timestep*femtosecond, {1:rpmd_beads[1]}) + else: + logger.info("Creating RPMD integrator with %i beads (no NB forces to contract).\n" % (rpmd_beads[0])) + integrator = RPMDIntegrator(rpmd_beads[0], temperature*kelvin, collision/picosecond, timestep*femtosecond) + else: + raise RuntimeError("Please provide a list of length 1, 2, or 3 to rpmd_beads") else: integrator = LangevinIntegrator(temperature*kelvin, collision/picosecond, timestep*femtosecond) else: ## If no temperature control, default to the Verlet integrator. - if rpmd_beads > 0: + if len(rpmd_beads) > 0: raise RuntimeError("No RPMD integrator without temperature control.") if mts: warn_once("No multiple timestep integrator without temperature control.") integrator = VerletIntegrator(timestep*femtoseconds) @@ -610,13 +644,13 @@ def create_simulation(self, timestep=1.0, faststep=0.25, temperature=None, press elif pressure != None: warn_once("Pressure is ignored because pbc is set to False.") ## Set up for energy component analysis. - if not mts: + if not mts and len(rpmd_beads) == 0: for i, j in enumerate(self.system.getForces()): j.setForceGroup(i) ## If virtual particles are used with AMOEBA... SetAmoebaVirtualExclusions(self.system) - + ## Finally create the simulation object. self.simulation = Simulation(self.mod.topology, self.system, integrator, self.platform) diff --git a/src/parser.py b/src/parser.py index 4cc510d4c..612ec3035 100644 --- a/src/parser.py +++ b/src/parser.py @@ -77,13 +77,13 @@ }, 'lists' : {"forcefield" : ([], 200, 'The names of force fields, corresponding to directory forcefields/file_name.(itp,xml,prm,frcmod,mol2)', 'All (important)'), "scanindex_num" : ([], -100, 'Numerical index of the parameter to scan over', 'Job types scan_mvals and scan_pvals'), - "scanindex_name" : ([], -100, 'Parameter name to scan over (should convert to a numerical index)', 'Job types scan_mvals and scan_pvals') + "scanindex_name" : ([], -100, 'Parameter name to scan over (should convert to a numerical index)', 'Job types scan_mvals and scan_pvals'), + "rpmd_beads" : ([], -160, 'Number of beads / NB contracted beads / Recip contracted beads in ring polymer MD', 'Condensed phase property targets (advanced usage)', 'liquid_openmm'), }, 'ints' : {"maxstep" : (100, 50, 'Maximum number of steps in an optimization', 'Main Optimizer'), "objective_history" : (2, 20, 'Number of good optimization steps to average over when checking the objective convergence criterion', 'Main Optimizer (jobtype "newton")'), "wq_port" : (0, 0, 'The port number to use for Work Queue', 'Targets that use Work Queue (advanced usage)'), "criteria" : (1, 160, 'The number of convergence criteria that must be met for main optimizer to converge', 'Main Optimizer'), - "rpmd_beads" : (0, -160, 'Number of beads in ring polymer MD (zero to disable)', 'Condensed phase property targets (advanced usage)', 'liquid_openmm'), }, 'bools' : {"backup" : (1, 10, 'Write temp directories to backup before wiping them'), "writechk_step" : (1, -50, 'Write the checkpoint file at every optimization step'), From 17ebcd3ab54fac01df23a26ffcbcf5220ced3760 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Wed, 19 Mar 2014 11:27:06 -0700 Subject: [PATCH 03/34] Added RPMD kinetic energy calculation --- src/openmmio.py | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/src/openmmio.py b/src/openmmio.py index 0d4129d0e..dc930707f 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -985,8 +985,24 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, if iteration >= 0: self.simulation.step(nsave) state = self.simulation.context.getState(getEnergy=True,getPositions=True,getVelocities=False,getForces=False) - kinetic = state.getKineticEnergy()/self.tdiv - potential = state.getPotentialEnergy() +##### + if self.tdiv>=2: + hbar=0.0635078*kilojoule*picosecond/mole + pimdstate=["?"]*self.tdiv + kinetic=0.0*kilojoule/mole + potential=0.0*kilojoule/mole + for i in range(self.tdiv): + pimdstate[i]=integrator.getstate(i,getEnergy=True,getPositions=True,group=-1) + potential=potential+pimdstate[i].getPotentialEnergy()/self.tdiv + kinetic=kinetic+pimdstate[i].getKineticEnergy()/self.tdiv + for i in range(self.tdiv): + ii=(i+1)%self.tdiv + deltabead=np.array(pimdstate[i].getPositions())-np.array(pimdstate[ii].getPositions()) + kinetic=kinetic-((deltabead*deltabead).sum(axis=1))*np.array(getParticleMass(j) for j in range system.getNumParticles())*(kB**2*integrator.getTemperature()**2*self.tdiv)/(2.0*hbar**2) + else: + kinetic = state.getKineticEnergy() + potential = state.getPotentialEnergy() +##### if self.pbc: box_vectors = state.getPeriodicBoxVectors() volume = self.compute_volume(box_vectors) From c5ce8b5b4a9b6d27330462dff007b059e209e892 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Wed, 19 Mar 2014 11:54:50 -0700 Subject: [PATCH 04/34] test commit --- src/openmmio.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/openmmio.py b/src/openmmio.py index dc930707f..9e852ea27 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -1002,7 +1002,7 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, else: kinetic = state.getKineticEnergy() potential = state.getPotentialEnergy() -##### +#####Energy Estimator if self.pbc: box_vectors = state.getPeriodicBoxVectors() volume = self.compute_volume(box_vectors) From 4186cabc9826b0471dae54f31e5b44f59d5a640f Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Wed, 19 Mar 2014 15:00:38 -0700 Subject: [PATCH 05/34] Test --- src/openmmio.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/openmmio.py b/src/openmmio.py index 9e852ea27..943e634f4 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -209,6 +209,9 @@ def CopyNonbondedParameters(src, dest): for i in range(src.getNumExceptions()): dest.setExceptionParameters(i,*src.getExceptionParameters(i)) +##def CopyCustomBondParameters(src, dest): +## + def do_nothing(src, dest): return @@ -226,6 +229,7 @@ def CopySystemParameters(src,dest): 'HarmonicAngleForce':CopyHarmonicAngleParameters, 'PeriodicTorsionForce':CopyPeriodicTorsionParameters, 'NonbondedForce':CopyNonbondedParameters, +## 'CustomBondForce':CopyCustomBondParameters, 'CMMotionRemover':do_nothing} for i in range(src.getNumForces()): nm = src.getForce(i).__class__.__name__ @@ -998,7 +1002,7 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, for i in range(self.tdiv): ii=(i+1)%self.tdiv deltabead=np.array(pimdstate[i].getPositions())-np.array(pimdstate[ii].getPositions()) - kinetic=kinetic-((deltabead*deltabead).sum(axis=1))*np.array(getParticleMass(j) for j in range system.getNumParticles())*(kB**2*integrator.getTemperature()**2*self.tdiv)/(2.0*hbar**2) + kinetic=kinetic-((deltabead*deltabead).sum(axis=1))*np.array([getParticleMass(j) for j in range system.getNumParticles()])*(kB**2*integrator.getTemperature()**2*self.tdiv)/(2.0*hbar**2) else: kinetic = state.getKineticEnergy() potential = state.getPotentialEnergy() From 98b2a156f0290d5d5183d299e92aa1a5a48e2cd6 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Wed, 19 Mar 2014 15:25:16 -0700 Subject: [PATCH 06/34] TEST --- src/openmmio.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/openmmio.py b/src/openmmio.py index 943e634f4..d9fedac39 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -1002,7 +1002,7 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, for i in range(self.tdiv): ii=(i+1)%self.tdiv deltabead=np.array(pimdstate[i].getPositions())-np.array(pimdstate[ii].getPositions()) - kinetic=kinetic-((deltabead*deltabead).sum(axis=1))*np.array([getParticleMass(j) for j in range system.getNumParticles()])*(kB**2*integrator.getTemperature()**2*self.tdiv)/(2.0*hbar**2) + kinetic=kinetic-((deltabead*deltabead).sum(axis=1))*np.array([getParticleMass(j) for j in range (system.getNumParticles())])*(kB**2*integrator.getTemperature()**2*self.tdiv)/(2.0*hbar**2) else: kinetic = state.getKineticEnergy() potential = state.getPotentialEnergy() From 7adaafeecfb3d4723e32f66ccdc74466d7efb8a5 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Wed, 19 Mar 2014 15:50:21 -0700 Subject: [PATCH 07/34] test1 --- src/openmmio.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/openmmio.py b/src/openmmio.py index d9fedac39..7c729e6c9 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -1002,7 +1002,7 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, for i in range(self.tdiv): ii=(i+1)%self.tdiv deltabead=np.array(pimdstate[i].getPositions())-np.array(pimdstate[ii].getPositions()) - kinetic=kinetic-((deltabead*deltabead).sum(axis=1))*np.array([getParticleMass(j) for j in range (system.getNumParticles())])*(kB**2*integrator.getTemperature()**2*self.tdiv)/(2.0*hbar**2) + kinetic=kinetic-np.sum(((deltabead*deltabead).sum(axis=1))*np.array([getParticleMass(j) for j in range (system.getNumParticles())]))*(kB**2*integrator.getTemperature()**2*self.tdiv)/(2.0*hbar**2) else: kinetic = state.getKineticEnergy() potential = state.getPotentialEnergy() From b0a74f983e199d59880cb59c26ba1c693cbac87c Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Wed, 19 Mar 2014 18:50:58 -0700 Subject: [PATCH 08/34] CustomBondForce --- src/openmmio.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/openmmio.py b/src/openmmio.py index 7c729e6c9..168d25182 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -209,9 +209,16 @@ def CopyNonbondedParameters(src, dest): for i in range(src.getNumExceptions()): dest.setExceptionParameters(i,*src.getExceptionParameters(i)) -##def CopyCustomBondParameters(src, dest): -## - +def CopyCustomBondParameters(src, dest): + dest.setEnergyFuncion(src.getEnergyFunction()) + for i in range(getNumGlobalParameters()): + dest.setGlobalParameterName(i,*src.GlobalParameterName(i)) + for i in range(getNumGlobalParameters()): + dest.setGlobalParameterDefaultValue(i,*src.getGlobalParameterDefaultValue(i)) + for i in range(getNumPerBondParameters()): + dest.setPerBondParameterName(i,*src.getPerBondParameterName(i)) + for i in range(getNumPerBondParameters()): + dest.setBondParameters(i,*src.getBondParameter(i)) def do_nothing(src, dest): return @@ -229,7 +236,7 @@ def CopySystemParameters(src,dest): 'HarmonicAngleForce':CopyHarmonicAngleParameters, 'PeriodicTorsionForce':CopyPeriodicTorsionParameters, 'NonbondedForce':CopyNonbondedParameters, -## 'CustomBondForce':CopyCustomBondParameters, + 'CustomBondForce':CopyCustomBondParameters, 'CMMotionRemover':do_nothing} for i in range(src.getNumForces()): nm = src.getForce(i).__class__.__name__ From ce5c26419f4c20e8e81fdf6e8dcc2241390f2767 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Thu, 20 Mar 2014 16:41:08 -0700 Subject: [PATCH 09/34] CustomBondForce(error in getPerBondParameterName) --- src/openmmio.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/openmmio.py b/src/openmmio.py index 168d25182..bf5b5f3bf 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -210,15 +210,15 @@ def CopyNonbondedParameters(src, dest): dest.setExceptionParameters(i,*src.getExceptionParameters(i)) def CopyCustomBondParameters(src, dest): - dest.setEnergyFuncion(src.getEnergyFunction()) - for i in range(getNumGlobalParameters()): + dest.setEnergyFunction(src.getEnergyFunction()) + for i in range(src.getNumGlobalParameters()): dest.setGlobalParameterName(i,*src.GlobalParameterName(i)) - for i in range(getNumGlobalParameters()): + for i in range(src.getNumGlobalParameters()): dest.setGlobalParameterDefaultValue(i,*src.getGlobalParameterDefaultValue(i)) - for i in range(getNumPerBondParameters()): + for i in range(src.getNumPerBondParameters()): dest.setPerBondParameterName(i,*src.getPerBondParameterName(i)) - for i in range(getNumPerBondParameters()): - dest.setBondParameters(i,*src.getBondParameter(i)) + for i in range(src.getNumPerBondParameters()): + dest.setBondParameters(i,*src.getBondParameters(i)) def do_nothing(src, dest): return From bee0c4dd5f01b4a17490e96106aebb3308088fe3 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Sat, 22 Mar 2014 16:45:35 -0700 Subject: [PATCH 10/34] CustomBondForce(working version) --- src/openmmio.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/openmmio.py b/src/openmmio.py index bf5b5f3bf..1b630b5d3 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -212,12 +212,12 @@ def CopyNonbondedParameters(src, dest): def CopyCustomBondParameters(src, dest): dest.setEnergyFunction(src.getEnergyFunction()) for i in range(src.getNumGlobalParameters()): - dest.setGlobalParameterName(i,*src.GlobalParameterName(i)) + dest.setGlobalParameterName(i,src.GlobalParameterName(i)) for i in range(src.getNumGlobalParameters()): dest.setGlobalParameterDefaultValue(i,*src.getGlobalParameterDefaultValue(i)) for i in range(src.getNumPerBondParameters()): - dest.setPerBondParameterName(i,*src.getPerBondParameterName(i)) - for i in range(src.getNumPerBondParameters()): + dest.setPerBondParameterName(i,src.getPerBondParameterName(i)) + for i in range(src.getNumBonds()): dest.setBondParameters(i,*src.getBondParameters(i)) def do_nothing(src, dest): return From d2a1c65e9d1bf7e23a984cb7578ff5b02ad4f556 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Sun, 23 Mar 2014 02:34:43 -0700 Subject: [PATCH 11/34] RPMD input should be list of integer not number --- src/openmmio.py | 59 ++++++++++++++++++++++++------------------------- 1 file changed, 29 insertions(+), 30 deletions(-) diff --git a/src/openmmio.py b/src/openmmio.py index 1b630b5d3..4dcbc32b0 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -584,7 +584,6 @@ def create_simulation(self, timestep=1.0, faststep=0.25, temperature=None, press options pertinent to system setup, including the type of MD integrator and type of pressure control. """ - # Divisor for the temperature (RPMD sets it to nonzero.) self.tdiv = 1 @@ -598,10 +597,10 @@ def create_simulation(self, timestep=1.0, faststep=0.25, temperature=None, press timestep*femtosecond, self.system, ninnersteps=int(timestep/faststep)) else: if len(rpmd_beads) > 0: - self.tdiv = rpmd_beads[0] + self.tdiv = int(rpmd_beads[0]) if len(rpmd_beads) == 1: - logger.info("Creating RPMD integrator with %i beads.\n" % rpmd_beads[0]) - integrator = RPMDIntegrator(rpmd_beads[0], temperature*kelvin, collision/picosecond, timestep*femtosecond) + logger.info("Creating RPMD integrator with %i beads.\n" % int(rpmd_beads[0])) + integrator = RPMDIntegrator(int(rpmd_beads[0]), temperature*kelvin, collision/picosecond, timestep*femtosecond) elif len(rpmd_beads) == 2: contract = False for frc in self.system.getForces(): @@ -609,11 +608,11 @@ def create_simulation(self, timestep=1.0, faststep=0.25, temperature=None, press contract = True frc.setForceGroup(1) if contract: - logger.info("Creating RPMD integrator with %i beads (NB forces contracted to %i).\n" % (rpmd_beads[0], rpmd_beads[1])) - integrator = RPMDIntegrator(rpmd_beads[0], temperature*kelvin, collision/picosecond, timestep*femtosecond, {1:rpmd_beads[1]}) + logger.info("Creating RPMD integrator with %i beads (NB forces contracted to %i).\n" % (int(rpmd_beads[0]), int(rpmd_beads[1]))) + integrator = RPMDIntegrator(int(rpmd_beads[0]), temperature*kelvin, collision/picosecond, timestep*femtosecond, {1:int(rpmd_beads[1])}) else: - logger.info("Creating RPMD integrator with %i beads (no NB forces to contract).\n" % (rpmd_beads[0])) - integrator = RPMDIntegrator(rpmd_beads[0], temperature*kelvin, collision/picosecond, timestep*femtosecond) + logger.info("Creating RPMD integrator with %i beads (no NB forces to contract).\n" % (int(rpmd_beads[0]))) + integrator = RPMDIntegrator(int(rpmd_beads[0]), temperature*kelvin, collision/picosecond, timestep*femtosecond) elif len(rpmd_beads) == 3: contract = False contract_recip = False @@ -625,14 +624,14 @@ def create_simulation(self, timestep=1.0, faststep=0.25, temperature=None, press contract_recip = True frc.setReciprocalSpaceForceGroup(2) if contract_recip: - logger.info("Creating RPMD integrator with %i beads (NB/Recip forces contracted to %i/%i).\n" % (rpmd_beads[0], rpmd_beads[1], rpmd_beads[2])) - integrator = RPMDIntegrator(rpmd_beads[0], temperature*kelvin, collision/picosecond, timestep*femtosecond, {1:rpmd_beads[1], 2:rpmd_beads[2]}) + logger.info("Creating RPMD integrator with %i beads (NB/Recip forces contracted to %i/%i).\n" % (int(rpmd_beads[0]), int(rpmd_beads[1]), int(rpmd_beads[2]))) + integrator = RPMDIntegrator(int(rpmd_beads[0]), temperature*kelvin, collision/picosecond, timestep*femtosecond, {1:int(rpmd_beads[1]), 2:int(rpmd_beads[2])}) elif contract: - logger.info("Creating RPMD integrator with %i beads (NB forces contracted to %i, no Recip).\n" % (rpmd_beads[0], rpmd_beads[1])) - integrator = RPMDIntegrator(rpmd_beads[0], temperature*kelvin, collision/picosecond, timestep*femtosecond, {1:rpmd_beads[1]}) + logger.info("Creating RPMD integrator with %i beads (NB forces contracted to %i, no Recip).\n" % (int(rpmd_beads[0]), int(rpmd_beads[1]))) + integrator = RPMDIntegrator(int(rpmd_beads[0]), temperature*kelvin, collision/picosecond, timestep*femtosecond, {1:int(rpmd_beads[1])}) else: - logger.info("Creating RPMD integrator with %i beads (no NB forces to contract).\n" % (rpmd_beads[0])) - integrator = RPMDIntegrator(rpmd_beads[0], temperature*kelvin, collision/picosecond, timestep*femtosecond) + logger.info("Creating RPMD integrator with %i beads (no NB forces to contract).\n" % (int(rpmd_beads[0]))) + integrator = RPMDIntegrator(int(rpmd_beads[0]), temperature*kelvin, collision/picosecond, timestep*femtosecond) else: raise RuntimeError("Please provide a list of length 1, 2, or 3 to rpmd_beads") else: @@ -997,22 +996,22 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, self.simulation.step(nsave) state = self.simulation.context.getState(getEnergy=True,getPositions=True,getVelocities=False,getForces=False) ##### - if self.tdiv>=2: - hbar=0.0635078*kilojoule*picosecond/mole - pimdstate=["?"]*self.tdiv - kinetic=0.0*kilojoule/mole - potential=0.0*kilojoule/mole - for i in range(self.tdiv): - pimdstate[i]=integrator.getstate(i,getEnergy=True,getPositions=True,group=-1) - potential=potential+pimdstate[i].getPotentialEnergy()/self.tdiv - kinetic=kinetic+pimdstate[i].getKineticEnergy()/self.tdiv - for i in range(self.tdiv): - ii=(i+1)%self.tdiv - deltabead=np.array(pimdstate[i].getPositions())-np.array(pimdstate[ii].getPositions()) - kinetic=kinetic-np.sum(((deltabead*deltabead).sum(axis=1))*np.array([getParticleMass(j) for j in range (system.getNumParticles())]))*(kB**2*integrator.getTemperature()**2*self.tdiv)/(2.0*hbar**2) - else: - kinetic = state.getKineticEnergy() - potential = state.getPotentialEnergy() +# if self.tdiv>=2: +# hbar=0.0635078*kilojoule*picosecond/mole +# pimdstate=["?"]*self.tdiv +# kinetic=0.0*kilojoule/mole +# potential=0.0*kilojoule/mole +# for i in range(self.tdiv): +# pimdstate[i]=integrator.getstate(i,getEnergy=True,getPositions=True,group=-1) +# potential=potential+pimdstate[i].getPotentialEnergy()/self.tdiv +# kinetic=kinetic+pimdstate[i].getKineticEnergy()/self.tdiv +# for i in range(self.tdiv): +# ii=(i+1)%self.tdiv +# deltabead=np.array(pimdstate[i].getPositions())-np.array(pimdstate[ii].getPositions()) +# kinetic=kinetic-np.sum(((deltabead*deltabead).sum(axis=1))*np.array([getParticleMass(j) for j in range (system.getNumParticles())]))*(kB**2*integrator.getTemperature()**2*self.tdiv)/(2.0*hbar**2) +# else: + kinetic = state.getKineticEnergy()/self.tdiv + potential = state.getPotentialEnergy() #####Energy Estimator if self.pbc: box_vectors = state.getPeriodicBoxVectors() From 4e15d7151f6cc5b54bb4ede1ef7b20c435e36814 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Wed, 26 Mar 2014 10:38:18 -0700 Subject: [PATCH 12/34] Energy estimator --- src/openmmio.py | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/src/openmmio.py b/src/openmmio.py index 4dcbc32b0..6623748c3 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -995,24 +995,24 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, if iteration >= 0: self.simulation.step(nsave) state = self.simulation.context.getState(getEnergy=True,getPositions=True,getVelocities=False,getForces=False) +#####Energy estimator + if self.tdiv>=2: + hbar=0.0635078*kilojoule*picosecond/mole + pimdstate=["?"]*self.tdiv + kinetic=0.0*kilojoule/mole + potential=0.0*kilojoule/mole + for i in range(self.tdiv): + pimdstate[i]=integrator.getstate(i,getEnergy=True,getPositions=True,group=-1) + potential=potential+pimdstate[i].getPotentialEnergy()/self.tdiv + kinetic=kinetic+pimdstate[i].getKineticEnergy()/self.tdiv + for i in range(self.tdiv): + ii=(i+1)%self.tdiv + deltabead=np.array(pimdstate[i].getPositions())-np.array(pimdstate[ii].getPositions()) + kinetic=kinetic-np.sum(((deltabead*deltabead).sum(axis=1))*np.array([getParticleMass(j) for j in range (system.getNumParticles())]))*(kB**2*integrator.getTemperature()**2*self.tdiv)/(2.0*hbar**2) + else: + kinetic = state.getKineticEnergy()/self.tdiv + potential = state.getPotentialEnergy() ##### -# if self.tdiv>=2: -# hbar=0.0635078*kilojoule*picosecond/mole -# pimdstate=["?"]*self.tdiv -# kinetic=0.0*kilojoule/mole -# potential=0.0*kilojoule/mole -# for i in range(self.tdiv): -# pimdstate[i]=integrator.getstate(i,getEnergy=True,getPositions=True,group=-1) -# potential=potential+pimdstate[i].getPotentialEnergy()/self.tdiv -# kinetic=kinetic+pimdstate[i].getKineticEnergy()/self.tdiv -# for i in range(self.tdiv): -# ii=(i+1)%self.tdiv -# deltabead=np.array(pimdstate[i].getPositions())-np.array(pimdstate[ii].getPositions()) -# kinetic=kinetic-np.sum(((deltabead*deltabead).sum(axis=1))*np.array([getParticleMass(j) for j in range (system.getNumParticles())]))*(kB**2*integrator.getTemperature()**2*self.tdiv)/(2.0*hbar**2) -# else: - kinetic = state.getKineticEnergy()/self.tdiv - potential = state.getPotentialEnergy() -#####Energy Estimator if self.pbc: box_vectors = state.getPeriodicBoxVectors() volume = self.compute_volume(box_vectors) From d2abf507aedd426716d8cd28c6d5f8b44950df2c Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Wed, 26 Mar 2014 14:15:22 -0700 Subject: [PATCH 13/34] Update rpmd energy estimating(centroid instead of primitive energy estimator) --- src/openmmio.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/openmmio.py b/src/openmmio.py index e7fea2b80..03927f6cb 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -1066,16 +1066,18 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, if self.tdiv>=2: hbar=0.0635078*kilojoule*picosecond/mole pimdstate=["?"]*self.tdiv + centroid=np.array([[0.0*nanometer,0.0*nanometer,0.0*nanometer]]*system.getNumParticles()) kinetic=0.0*kilojoule/mole potential=0.0*kilojoule/mole for i in range(self.tdiv): pimdstate[i]=integrator.getstate(i,getEnergy=True,getPositions=True,group=-1) + centroid=centroid+np.array(pimdstate[i].getPositions())/self.tdiv potential=potential+pimdstate[i].getPotentialEnergy()/self.tdiv - kinetic=kinetic+pimdstate[i].getKineticEnergy()/self.tdiv + kinetic=kinetic+pimdstate[i].getKineticEnergy()/(self.tdiv*self.tdiv) for i in range(self.tdiv): - ii=(i+1)%self.tdiv - deltabead=np.array(pimdstate[i].getPositions())-np.array(pimdstate[ii].getPositions()) - kinetic=kinetic-np.sum(((deltabead*deltabead).sum(axis=1))*np.array([getParticleMass(j) for j in range (system.getNumParticles())]))*(kB**2*integrator.getTemperature()**2*self.tdiv)/(2.0*hbar**2) + dif=np.array(pimdstate[i].getPositions())-centroid + der=-1.0*np.array(pimdstates[i].getForces()) + kinetic=kinetic+np.sum((dif*der).sum(axis=1))*0.5/self.tdiv else: kinetic = state.getKineticEnergy()/self.tdiv potential = state.getPotentialEnergy() From 16939e4a853d5c5b0daabddf744cd7d42cbdcf96 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Thu, 27 Mar 2014 13:39:36 -0700 Subject: [PATCH 14/34] Centroid Energy Estimator --- src/openmmio.py | 41 +++++++++++++++++++++++++++-------------- 1 file changed, 27 insertions(+), 14 deletions(-) diff --git a/src/openmmio.py b/src/openmmio.py index 03927f6cb..1b3e07eaf 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -696,7 +696,6 @@ def create_simulation(self, timestep=1.0, faststep=0.25, temperature=None, press raise RuntimeError("No RPMD integrator without temperature control.") if mts: warn_once("No multiple timestep integrator without temperature control.") integrator = VerletIntegrator(timestep*femtoseconds) - ## Add the barostat. if pressure != None: if anisotropic: @@ -717,13 +716,11 @@ def create_simulation(self, timestep=1.0, faststep=0.25, temperature=None, press ## Finally create the simulation object. self.simulation = Simulation(self.mod.topology, self.system, integrator, self.platform) - ## Print platform properties. # logger.info("I'm using the platform %s\n" % self.simulation.context.getPlatform().getName()) # printcool_dictionary({i:self.simulation.context.getPlatform().getPropertyValue(self.simulation.context,i) \ # for i in self.simulation.context.getPlatform().getPropertyNames()}, \ # title="Platform %s has properties:" % self.simulation.context.getPlatform().getName()) - def update_simulation(self, **kwargs): """ @@ -758,7 +755,6 @@ def update_simulation(self, **kwargs): UpdateSimulationParameters(self.system, self.simulation) else: self.create_simulation(**self.simkwargs) - def set_positions(self, shot=0, traj=None): """ @@ -784,7 +780,6 @@ def set_positions(self, shot=0, traj=None): # self.simulation.context.setPositions(ResetVirtualSites_fast(self.xyz_omms[shot][0], self.vsinfo)) self.simulation.context.setPositions(self.xyz_omms[shot][0]) self.simulation.context.computeVirtualSites() - def compute_volume(self, box_vectors): """ Compute the total volume of an OpenMM system. """ [a,b,c] = box_vectors @@ -972,7 +967,6 @@ def interaction_energy(self, fraga, fragb): B = self.B.energy() return (D - A - B) / 4.184 - def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, nequil=0, nsave=1000, minimize=True, anisotropic=False, save_traj=False, verbose=False, **kwargs): """ @@ -1034,7 +1028,6 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, # Determine number of degrees of freedom. self.ndof = 3*(self.system.getNumParticles() - sum([self.system.isVirtualSite(i) for i in range(self.system.getNumParticles())])) \ - self.system.getNumConstraints() - 3*self.pbc - # Initialize statistics. edecomp = OrderedDict() # Stored coordinates, box vectors @@ -1065,21 +1058,22 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, #####Energy estimator if self.tdiv>=2: hbar=0.0635078*kilojoule*picosecond/mole - pimdstate=["?"]*self.tdiv - centroid=np.array([[0.0*nanometer,0.0*nanometer,0.0*nanometer]]*system.getNumParticles()) + pimdstate=[] + centroid=np.array([[0.0*nanometer,0.0*nanometer,0.0*nanometer]]*self.system.getNumParticles()) kinetic=0.0*kilojoule/mole potential=0.0*kilojoule/mole + integrator=self.simulation.context.getIntegrator() for i in range(self.tdiv): - pimdstate[i]=integrator.getstate(i,getEnergy=True,getPositions=True,group=-1) + pimdstate.append(integrator.getState(i,getEnergy=True,getPositions=True,getForces=True,groups=-1)) centroid=centroid+np.array(pimdstate[i].getPositions())/self.tdiv potential=potential+pimdstate[i].getPotentialEnergy()/self.tdiv kinetic=kinetic+pimdstate[i].getKineticEnergy()/(self.tdiv*self.tdiv) for i in range(self.tdiv): dif=np.array(pimdstate[i].getPositions())-centroid - der=-1.0*np.array(pimdstates[i].getForces()) + der=-1.0*np.array(pimdstate[i].getForces()) kinetic=kinetic+np.sum((dif*der).sum(axis=1))*0.5/self.tdiv else: - kinetic = state.getKineticEnergy()/self.tdiv + kinetic = state.getKineticEnergy() potential = state.getPotentialEnergy() ##### if self.pbc: @@ -1111,8 +1105,27 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, if iteration >= 0: self.simulation.step(nsave) # Compute properties. state = self.simulation.context.getState(getEnergy=True,getPositions=True,getVelocities=False,getForces=False) - kinetic = state.getKineticEnergy()/self.tdiv - potential = state.getPotentialEnergy() +#####RPMD Energy estimator + if self.tdiv>=2: + hbar=0.0635078*kilojoule*picosecond/mole + pimdstate=[] + centroid=np.array([[0.0*nanometer,0.0*nanometer,0.0*nanometer]]*self.system.getNumParticles()) + kinetic=0.0*kilojoule/mole + potential=0.0*kilojoule/mole + integrator=self.simulation.context.getIntegrator() + for i in range(self.tdiv): + pimdstate.append(integrator.getState(i,getEnergy=True,getPositions=True,getForces=True,groups=-1)) + centroid=centroid+np.array(pimdstate[i].getPositions())/self.tdiv + potential=potential+pimdstate[i].getPotentialEnergy()/self.tdiv + kinetic=kinetic+pimdstate[i].getKineticEnergy()/(self.tdiv*self.tdiv) + for i in range(self.tdiv): + dif=np.array(pimdstate[i].getPositions())-centroid + der=-1.0*np.array(pimdstate[i].getForces()) + kinetic=kinetic+np.sum((dif*der).sum(axis=1))*0.5/self.tdiv + else: + kinetic = state.getKineticEnergy() + potential = state.getPotentialEnergy() +##### kinetic_temperature = 2.0 * kinetic / kB / self.ndof if self.pbc: box_vectors = state.getPeriodicBoxVectors() From 0ed58f53d88b539b2804a5eeb94713c814c36b04 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Wed, 2 Apr 2014 13:27:17 -0700 Subject: [PATCH 15/34] test setposition --- src/openmmio.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/openmmio.py b/src/openmmio.py index 1b3e07eaf..e8938ae00 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -774,12 +774,19 @@ def set_positions(self, shot=0, traj=None): # simulation.context.computeVirtualSites() #---- # NOTE: Periodic box vectors must be set FIRST - if self.pbc: - self.simulation.context.setPeriodicBoxVectors(*self.xyz_omms[shot][1]) + if self.tdiv==1: + if self.pbc: + self.simulation.context.setPeriodicBoxVectors(*self.xyz_omms[shot][1]) # self.simulation.context.setPositions(ResetVirtualSites(self.xyz_omms[shot][0], self.system)) # self.simulation.context.setPositions(ResetVirtualSites_fast(self.xyz_omms[shot][0], self.vsinfo)) - self.simulation.context.setPositions(self.xyz_omms[shot][0]) - self.simulation.context.computeVirtualSites() + self.simulation.context.setPositions(self.xyz_omms[shot][0]) + self.simulation.context.computeVirtualSites() + if self.tdiv>1: + if self.pbc: + self.simulation.context.setPeriodicBoxVectors(*self.xyz_omms[shot][1]) + for i in range(self.tdiv): + self.simulation.context.getIntegrator.setPositions(i,self.simulation.context.getIngegrator.getState(i,getPositions=True).getPositions()) + def compute_volume(self, box_vectors): """ Compute the total volume of an OpenMM system. """ [a,b,c] = box_vectors @@ -1057,7 +1064,6 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, state = self.simulation.context.getState(getEnergy=True,getPositions=True,getVelocities=False,getForces=False) #####Energy estimator if self.tdiv>=2: - hbar=0.0635078*kilojoule*picosecond/mole pimdstate=[] centroid=np.array([[0.0*nanometer,0.0*nanometer,0.0*nanometer]]*self.system.getNumParticles()) kinetic=0.0*kilojoule/mole From 1c555829af3e7ed60c15f07d7c6bd801922654ea Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Fri, 4 Apr 2014 15:58:42 -0700 Subject: [PATCH 16/34] RPMD set_position (test) --- src/openmmio.py | 156 +++++++++++++++++++++++++++++++----------------- 1 file changed, 102 insertions(+), 54 deletions(-) diff --git a/src/openmmio.py b/src/openmmio.py index e8938ae00..61dd953fb 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -54,50 +54,96 @@ def get_multipoles(simulation,q=None,mass=None,positions=None,rmcom=True): qyz = 0.0 qzz = 0.0 enm_debye = 48.03204255928332 # Conversion factor from e*nm to Debye - for i in simulation.system.getForces(): - if isinstance(i, AmoebaMultipoleForce): - mm = i.getSystemMultipoleMoments(simulation.context) - dx += mm[1] - dy += mm[2] - dz += mm[3] - qxx += mm[4] - qxy += mm[5] - qxz += mm[6] - qyy += mm[8] - qyz += mm[9] - qzz += mm[12] - if isinstance(i, NonbondedForce): - # Get array of charges. - if q == None: - q = np.array([i.getParticleParameters(j)[0]._value for j in range(i.getNumParticles())]) - # Get array of positions in nanometers. - if positions == None: - positions = simulation.context.getState(getPositions=True).getPositions() - if mass == None: - mass = np.array([simulation.context.getSystem().getParticleMass(k).value_in_unit(dalton) \ - for k in range(simulation.context.getSystem().getNumParticles())]) - x = np.array(positions.value_in_unit(nanometer)) - if rmcom: - com = np.sum(x*mass.reshape(-1,1),axis=0) / np.sum(mass) - x -= com - xx, xy, xz, yy, yz, zz = (x[:,k]*x[:,l] for k, l in [(0,0),(0,1),(0,2),(1,1),(1,2),(2,2)]) - # Multiply charges by positions to get dipole moment. - dip = enm_debye * np.sum(x*q.reshape(-1,1),axis=0) - dx += dip[0] - dy += dip[1] - dz += dip[2] - qxx += enm_debye * 15 * np.sum(q*xx) - qxy += enm_debye * 15 * np.sum(q*xy) - qxz += enm_debye * 15 * np.sum(q*xz) - qyy += enm_debye * 15 * np.sum(q*yy) - qyz += enm_debye * 15 * np.sum(q*yz) - qzz += enm_debye * 15 * np.sum(q*zz) - tr = qxx+qyy+qzz - qxx -= tr/3 - qyy -= tr/3 - qzz -= tr/3 - # This ordering has to do with the way TINKER prints it out. - return [dx,dy,dz,qxx,qxy,qyy,qxz,qyz,qzz] + if(str(simulation.context.getIntegrator())[21:25]!='RPMD'):#if not RPMD integrator^-^ + for i in simulation.system.getForces(): + if isinstance(i, AmoebaMultipoleForce): + mm = i.getSystemMultipoleMoments(simulation.context) + dx += mm[1] + dy += mm[2] + dz += mm[3] + qxx += mm[4] + qxy += mm[5] + qxz += mm[6] + qyy += mm[8] + qyz += mm[9] + qzz += mm[12] + if isinstance(i, NonbondedForce): + # Get array of charges. + if q == None: + q = np.array([i.getParticleParameters(j)[0]._value for j in range(i.getNumParticles())]) + # Get array of positions in nanometers. + if positions == None: + positions = simulation.context.getState(getPositions=True).getPositions() + if mass == None: + mass = np.array([simulation.context.getSystem().getParticleMass(k).value_in_unit(dalton) \ + for k in range(simulation.context.getSystem().getNumParticles())]) + x = np.array(positions.value_in_unit(nanometer)) + if rmcom: + com = np.sum(x*mass.reshape(-1,1),axis=0) / np.sum(mass) + x -= com + xx, xy, xz, yy, yz, zz = (x[:,k]*x[:,l] for k, l in [(0,0),(0,1),(0,2),(1,1),(1,2),(2,2)]) + # Multiply charges by positions to get dipole moment. + dip = enm_debye * np.sum(x*q.reshape(-1,1),axis=0) + dx += dip[0] + dy += dip[1] + dz += dip[2] + qxx += enm_debye * 15 * np.sum(q*xx) + qxy += enm_debye * 15 * np.sum(q*xy) + qxz += enm_debye * 15 * np.sum(q*xz) + qyy += enm_debye * 15 * np.sum(q*yy) + qyz += enm_debye * 15 * np.sum(q*yz) + qzz += enm_debye * 15 * np.sum(q*zz) + tr = qxx+qyy+qzz + qxx -= tr/3 + qyy -= tr/3 + qzz -= tr/3 + # This ordering has to do with the way TINKER prints it out. + return [dx,dy,dz,qxx,qxy,qyy,qxz,qyz,qzz] + else:#RPMD integrator^-^(no edit yet) + for i in simulation.system.getForces(): + if isinstance(i, AmoebaMultipoleForce): + mm = i.getSystemMultipoleMoments(simulation.context) + dx += mm[1] + dy += mm[2] + dz += mm[3] + qxx += mm[4] + qxy += mm[5] + qxz += mm[6] + qyy += mm[8] + qyz += mm[9] + qzz += mm[12] + if isinstance(i, NonbondedForce): + # Get array of charges. + if q == None: + q = np.array([i.getParticleParameters(j)[0]._value for j in range(i.getNumParticles())]) + # Get array of positions in nanometers. + if positions == None: + positions = simulation.context.getState(getPositions=True).getPositions() + if mass == None: + mass = np.array([simulation.context.getSystem().getParticleMass(k).value_in_unit(dalton) \ + for k in range(simulation.context.getSystem().getNumParticles())]) + x = np.array(positions.value_in_unit(nanometer)) + if rmcom: + com = np.sum(x*mass.reshape(-1,1),axis=0) / np.sum(mass) + x -= com + xx, xy, xz, yy, yz, zz = (x[:,k]*x[:,l] for k, l in [(0,0),(0,1),(0,2),(1,1),(1,2),(2,2)]) + # Multiply charges by positions to get dipole moment. + dip = enm_debye * np.sum(x*q.reshape(-1,1),axis=0) + dx += dip[0] + dy += dip[1] + dz += dip[2] + qxx += enm_debye * 15 * np.sum(q*xx) + qxy += enm_debye * 15 * np.sum(q*xy) + qxz += enm_debye * 15 * np.sum(q*xz) + qyy += enm_debye * 15 * np.sum(q*yy) + qyz += enm_debye * 15 * np.sum(q*yz) + qzz += enm_debye * 15 * np.sum(q*zz) + tr = qxx+qyy+qzz + qxx -= tr/3 + qyy -= tr/3 + qzz -= tr/3 + # This ordering has to do with the way TINKER prints it out. + return [dx,dy,dz,qxx,qxy,qyy,qxz,qyz,qzz] def get_dipole(simulation,q=None,mass=None,positions=None): """Return the current dipole moment in Debye. @@ -774,18 +820,17 @@ def set_positions(self, shot=0, traj=None): # simulation.context.computeVirtualSites() #---- # NOTE: Periodic box vectors must be set FIRST - if self.tdiv==1: - if self.pbc: - self.simulation.context.setPeriodicBoxVectors(*self.xyz_omms[shot][1]) + if self.pbc: + self.simulation.context.setPeriodicBoxVectors(*self.xyz_omms[shot][1]) # self.simulation.context.setPositions(ResetVirtualSites(self.xyz_omms[shot][0], self.system)) # self.simulation.context.setPositions(ResetVirtualSites_fast(self.xyz_omms[shot][0], self.vsinfo)) - self.simulation.context.setPositions(self.xyz_omms[shot][0]) - self.simulation.context.computeVirtualSites() + self.simulation.context.setPositions(self.xyz_omms[shot][0]) + self.simulation.context.computeVirtualSites() if self.tdiv>1: - if self.pbc: - self.simulation.context.setPeriodicBoxVectors(*self.xyz_omms[shot][1]) for i in range(self.tdiv): - self.simulation.context.getIntegrator.setPositions(i,self.simulation.context.getIngegrator.getState(i,getPositions=True).getPositions()) + integrator=self.simulation.context.getIntegrator() + integrator.setPositions(i,self.xyz_omms[shot][0]) + ##compute virtual site method?^-^ def compute_volume(self, box_vectors): """ Compute the total volume of an OpenMM system. """ @@ -1051,6 +1096,7 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, #========================# # Initialize velocities. self.simulation.context.setVelocitiesToTemperature(temperature*kelvin) + # Equilibrate. if iequil > 0: if verbose: logger.info("Equilibrating...\n") @@ -1071,6 +1117,7 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, integrator=self.simulation.context.getIntegrator() for i in range(self.tdiv): pimdstate.append(integrator.getState(i,getEnergy=True,getPositions=True,getForces=True,groups=-1)) + for i in range(self.tdiv): centroid=centroid+np.array(pimdstate[i].getPositions())/self.tdiv potential=potential+pimdstate[i].getPotentialEnergy()/self.tdiv kinetic=kinetic+pimdstate[i].getKineticEnergy()/(self.tdiv*self.tdiv) @@ -1079,7 +1126,7 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, der=-1.0*np.array(pimdstate[i].getForces()) kinetic=kinetic+np.sum((dif*der).sum(axis=1))*0.5/self.tdiv else: - kinetic = state.getKineticEnergy() + kinetic = state.getKineticEnergy()/self.tdiv potential = state.getPotentialEnergy() ##### if self.pbc: @@ -1121,6 +1168,7 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, integrator=self.simulation.context.getIntegrator() for i in range(self.tdiv): pimdstate.append(integrator.getState(i,getEnergy=True,getPositions=True,getForces=True,groups=-1)) + for i in range(self.tdiv): centroid=centroid+np.array(pimdstate[i].getPositions())/self.tdiv potential=potential+pimdstate[i].getPotentialEnergy()/self.tdiv kinetic=kinetic+pimdstate[i].getKineticEnergy()/(self.tdiv*self.tdiv) @@ -1129,7 +1177,7 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, der=-1.0*np.array(pimdstate[i].getForces()) kinetic=kinetic+np.sum((dif*der).sum(axis=1))*0.5/self.tdiv else: - kinetic = state.getKineticEnergy() + kinetic = state.getKineticEnergy()/self.tdiv potential = state.getPotentialEnergy() ##### kinetic_temperature = 2.0 * kinetic / kB / self.ndof From c8524e41857d8cf2cc6f010e250497322b5957b0 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Fri, 4 Apr 2014 17:25:46 -0700 Subject: [PATCH 17/34] Undo get multipole function --- src/openmmio.py | 134 ++++++++++++++++-------------------------------- 1 file changed, 44 insertions(+), 90 deletions(-) diff --git a/src/openmmio.py b/src/openmmio.py index 61dd953fb..78c00bf37 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -54,96 +54,50 @@ def get_multipoles(simulation,q=None,mass=None,positions=None,rmcom=True): qyz = 0.0 qzz = 0.0 enm_debye = 48.03204255928332 # Conversion factor from e*nm to Debye - if(str(simulation.context.getIntegrator())[21:25]!='RPMD'):#if not RPMD integrator^-^ - for i in simulation.system.getForces(): - if isinstance(i, AmoebaMultipoleForce): - mm = i.getSystemMultipoleMoments(simulation.context) - dx += mm[1] - dy += mm[2] - dz += mm[3] - qxx += mm[4] - qxy += mm[5] - qxz += mm[6] - qyy += mm[8] - qyz += mm[9] - qzz += mm[12] - if isinstance(i, NonbondedForce): - # Get array of charges. - if q == None: - q = np.array([i.getParticleParameters(j)[0]._value for j in range(i.getNumParticles())]) - # Get array of positions in nanometers. - if positions == None: - positions = simulation.context.getState(getPositions=True).getPositions() - if mass == None: - mass = np.array([simulation.context.getSystem().getParticleMass(k).value_in_unit(dalton) \ - for k in range(simulation.context.getSystem().getNumParticles())]) - x = np.array(positions.value_in_unit(nanometer)) - if rmcom: - com = np.sum(x*mass.reshape(-1,1),axis=0) / np.sum(mass) - x -= com - xx, xy, xz, yy, yz, zz = (x[:,k]*x[:,l] for k, l in [(0,0),(0,1),(0,2),(1,1),(1,2),(2,2)]) - # Multiply charges by positions to get dipole moment. - dip = enm_debye * np.sum(x*q.reshape(-1,1),axis=0) - dx += dip[0] - dy += dip[1] - dz += dip[2] - qxx += enm_debye * 15 * np.sum(q*xx) - qxy += enm_debye * 15 * np.sum(q*xy) - qxz += enm_debye * 15 * np.sum(q*xz) - qyy += enm_debye * 15 * np.sum(q*yy) - qyz += enm_debye * 15 * np.sum(q*yz) - qzz += enm_debye * 15 * np.sum(q*zz) - tr = qxx+qyy+qzz - qxx -= tr/3 - qyy -= tr/3 - qzz -= tr/3 - # This ordering has to do with the way TINKER prints it out. - return [dx,dy,dz,qxx,qxy,qyy,qxz,qyz,qzz] - else:#RPMD integrator^-^(no edit yet) - for i in simulation.system.getForces(): - if isinstance(i, AmoebaMultipoleForce): - mm = i.getSystemMultipoleMoments(simulation.context) - dx += mm[1] - dy += mm[2] - dz += mm[3] - qxx += mm[4] - qxy += mm[5] - qxz += mm[6] - qyy += mm[8] - qyz += mm[9] - qzz += mm[12] - if isinstance(i, NonbondedForce): - # Get array of charges. - if q == None: - q = np.array([i.getParticleParameters(j)[0]._value for j in range(i.getNumParticles())]) - # Get array of positions in nanometers. - if positions == None: - positions = simulation.context.getState(getPositions=True).getPositions() - if mass == None: - mass = np.array([simulation.context.getSystem().getParticleMass(k).value_in_unit(dalton) \ - for k in range(simulation.context.getSystem().getNumParticles())]) - x = np.array(positions.value_in_unit(nanometer)) - if rmcom: - com = np.sum(x*mass.reshape(-1,1),axis=0) / np.sum(mass) - x -= com - xx, xy, xz, yy, yz, zz = (x[:,k]*x[:,l] for k, l in [(0,0),(0,1),(0,2),(1,1),(1,2),(2,2)]) - # Multiply charges by positions to get dipole moment. - dip = enm_debye * np.sum(x*q.reshape(-1,1),axis=0) - dx += dip[0] - dy += dip[1] - dz += dip[2] - qxx += enm_debye * 15 * np.sum(q*xx) - qxy += enm_debye * 15 * np.sum(q*xy) - qxz += enm_debye * 15 * np.sum(q*xz) - qyy += enm_debye * 15 * np.sum(q*yy) - qyz += enm_debye * 15 * np.sum(q*yz) - qzz += enm_debye * 15 * np.sum(q*zz) - tr = qxx+qyy+qzz - qxx -= tr/3 - qyy -= tr/3 - qzz -= tr/3 - # This ordering has to do with the way TINKER prints it out. - return [dx,dy,dz,qxx,qxy,qyy,qxz,qyz,qzz] + for i in simulation.system.getForces(): + if isinstance(i, AmoebaMultipoleForce): + mm = i.getSystemMultipoleMoments(simulation.context) + dx += mm[1] + dy += mm[2] + dz += mm[3] + qxx += mm[4] + qxy += mm[5] + qxz += mm[6] + qyy += mm[8] + qyz += mm[9] + qzz += mm[12] + if isinstance(i, NonbondedForce): + # Get array of charges. + if q == None: + q = np.array([i.getParticleParameters(j)[0]._value for j in range(i.getNumParticles())]) + # Get array of positions in nanometers. + if positions == None: + positions = simulation.context.getState(getPositions=True).getPositions() + if mass == None: + mass = np.array([simulation.context.getSystem().getParticleMass(k).value_in_unit(dalton) \ + for k in range(simulation.context.getSystem().getNumParticles())]) + x = np.array(positions.value_in_unit(nanometer)) + if rmcom: + com = np.sum(x*mass.reshape(-1,1),axis=0) / np.sum(mass) + x -= com + xx, xy, xz, yy, yz, zz = (x[:,k]*x[:,l] for k, l in [(0,0),(0,1),(0,2),(1,1),(1,2),(2,2)]) + # Multiply charges by positions to get dipole moment. + dip = enm_debye * np.sum(x*q.reshape(-1,1),axis=0) + dx += dip[0] + dy += dip[1] + dz += dip[2] + qxx += enm_debye * 15 * np.sum(q*xx) + qxy += enm_debye * 15 * np.sum(q*xy) + qxz += enm_debye * 15 * np.sum(q*xz) + qyy += enm_debye * 15 * np.sum(q*yy) + qyz += enm_debye * 15 * np.sum(q*yz) + qzz += enm_debye * 15 * np.sum(q*zz) + tr = qxx+qyy+qzz + qxx -= tr/3 + qyy -= tr/3 + qzz -= tr/3 + # This ordering has to do with the way TINKER prints it out. + return [dx,dy,dz,qxx,qxy,qyy,qxz,qyz,qzz] def get_dipole(simulation,q=None,mass=None,positions=None): """Return the current dipole moment in Debye. From 5ffd3fb5155759b331011a83ec6d62cc2c9475f2 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Mon, 7 Apr 2014 12:32:37 -0700 Subject: [PATCH 18/34] update energy estimator --- src/openmmio.py | 52 ++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 41 insertions(+), 11 deletions(-) diff --git a/src/openmmio.py b/src/openmmio.py index 78c00bf37..603a222a9 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -1062,26 +1062,42 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, if iteration >= 0: self.simulation.step(nsave) state = self.simulation.context.getState(getEnergy=True,getPositions=True,getVelocities=False,getForces=False) -#####Energy estimator + +##### Virial Energy estimator if self.tdiv>=2: + hbar=0.0635078*nanometer**2*dalton/picosecond pimdstate=[] centroid=np.array([[0.0*nanometer,0.0*nanometer,0.0*nanometer]]*self.system.getNumParticles()) kinetic=0.0*kilojoule/mole + cent_kinetic=0.0*kilojoule/mole + pri_kinetic=0.0*kilojoule/mole potential=0.0*kilojoule/mole integrator=self.simulation.context.getIntegrator() for i in range(self.tdiv): pimdstate.append(integrator.getState(i,getEnergy=True,getPositions=True,getForces=True,groups=-1)) for i in range(self.tdiv): centroid=centroid+np.array(pimdstate[i].getPositions())/self.tdiv - potential=potential+pimdstate[i].getPotentialEnergy()/self.tdiv - kinetic=kinetic+pimdstate[i].getKineticEnergy()/(self.tdiv*self.tdiv) - for i in range(self.tdiv): + potential=potential+pimdstate[i].getPotentialEnergy()/self.tdiv#potential energy + cent_kinetic=cent_kinetic+pimdstate[i].getKineticEnergy()/(self.tdiv*self.tdiv)#1st term in centroid energy estimator + pri_kinetic=pri_kinetic+pimdstate[i].getKineticEnergy()/self.tdiv#1st term in primitive kinetic energy estimator + kinetic=kinetic+pimdstate[i].getKineticEnergy()/(self.tdiv*self.tdiv)#the kinetic energy that average to thermostat temperature + for i in range(self.tdiv):#2nd term in centroid kinetic energy estimator dif=np.array(pimdstate[i].getPositions())-centroid der=-1.0*np.array(pimdstate[i].getForces()) - kinetic=kinetic+np.sum((dif*der).sum(axis=1))*0.5/self.tdiv + cent_kinetic=cent_kinetic+np.sum((dif*der).sum(axis=1))*0.5/self.tdiv + mass_matrix=[] + for i in range(self.system.getNumParticles()): + mass_matrix.append(self.system.getParticleMass(i)) + mass_matrix=np.array(mass_matrix) + for i in range(self.tdiv):#2nd term in primitive kinetic energy estimator + j=(i+1)%(self.tdiv) + beaddif=np.array(pimdstate[j].getPositions())-np.array(pimdstate[i].getPositions())#difference in position between ith and i+1th bead + pri_kinetic=pri_kinetic-np.sum(((beaddif*beaddif).sum(axis=1))*mass_matrix*(kB**2*integrator.getTemperature()**2*self.tdiv)/(2.0*hbar**2)) else: kinetic = state.getKineticEnergy()/self.tdiv potential = state.getPotentialEnergy() + cent_kinetic=0.0*kilojoule/mole + pri_kinetic=0.0*kilojoule/mole#Don't use centroid or primitive kinetic energy estimator for non-RPMD simulation ##### if self.pbc: box_vectors = state.getPeriodicBoxVectors() @@ -1112,27 +1128,41 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, if iteration >= 0: self.simulation.step(nsave) # Compute properties. state = self.simulation.context.getState(getEnergy=True,getPositions=True,getVelocities=False,getForces=False) -#####RPMD Energy estimator +#####Virial Energy estimator if self.tdiv>=2: - hbar=0.0635078*kilojoule*picosecond/mole + hbar=0.0635078*nanometer**2*dalton/picosecond pimdstate=[] centroid=np.array([[0.0*nanometer,0.0*nanometer,0.0*nanometer]]*self.system.getNumParticles()) kinetic=0.0*kilojoule/mole + cent_kinetic=0.0*kilojoule/mole + pri_kinetic=0.0*kilojoule/mole potential=0.0*kilojoule/mole integrator=self.simulation.context.getIntegrator() for i in range(self.tdiv): pimdstate.append(integrator.getState(i,getEnergy=True,getPositions=True,getForces=True,groups=-1)) for i in range(self.tdiv): centroid=centroid+np.array(pimdstate[i].getPositions())/self.tdiv - potential=potential+pimdstate[i].getPotentialEnergy()/self.tdiv - kinetic=kinetic+pimdstate[i].getKineticEnergy()/(self.tdiv*self.tdiv) - for i in range(self.tdiv): + potential=potential+pimdstate[i].getPotentialEnergy()/self.tdiv#potential energy + cent_kinetic=cent_kinetic+pimdstate[i].getKineticEnergy()/(self.tdiv*self.tdiv)#1st term in centroid energy estimator + pri_kinetic=pri_kinetic+pimdstate[i].getKineticEnergy()/self.tdiv#1st term in primitive kinetic energy estimator + kinetic=kinetic+pimdstate[i].getKineticEnergy()/(self.tdiv*self.tdiv)#the kinetic energy that average to thermostat temperature + for i in range(self.tdiv):#2nd term in centroid kinetic energy estimator dif=np.array(pimdstate[i].getPositions())-centroid der=-1.0*np.array(pimdstate[i].getForces()) - kinetic=kinetic+np.sum((dif*der).sum(axis=1))*0.5/self.tdiv + cent_kinetic=cent_kinetic+np.sum((dif*der).sum(axis=1))*0.5/self.tdiv + mass_matrix=[] + for i in range(self.system.getNumParticles()): + mass_matrix.append(self.system.getParticleMass(i)) + mass_matrix=np.array(mass_matrix) + for i in range(self.tdiv):#2nd term in primitive kinetic energy estimator + j=(i+1)%(self.tdiv) + beaddif=np.array(pimdstate[j].getPositions())-np.array(pimdstate[i].getPositions())#difference in position between ith and i+1th bead + pri_kinetic=pri_kinetic-np.sum(((beaddif*beaddif).sum(axis=1))*mass_matrix*(kB**2*integrator.getTemperature()**2*self.tdiv)/(2.0*hbar**2)) else: kinetic = state.getKineticEnergy()/self.tdiv potential = state.getPotentialEnergy() + cent_kinetic=0.0*kilojoule/mole + pri_kinetic=0.0*kilojoule/mole#Don't use centroid or primitive kinetic energy estimator for non-RPMD simulation ##### kinetic_temperature = 2.0 * kinetic / kB / self.ndof if self.pbc: From 5f55c7e0ccaa211d52ca02d8e970697a716b7f0f Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Tue, 8 Apr 2014 17:13:04 -0700 Subject: [PATCH 19/34] Implement Potential and Kinetic energy calculators for RPMD simulation as function --- src/openmmio.py | 53 ++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 52 insertions(+), 1 deletion(-) diff --git a/src/openmmio.py b/src/openmmio.py index 603a222a9..228fa6e44 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -41,7 +41,58 @@ def energy_components(Sim, verbose=False): for i in range(Sim.system.getNumForces()): EnergyTerms[Sim.system.getForce(i).__class__.__name__] = Sim.context.getState(getEnergy=True,groups=2**i).getPotentialEnergy() / kilojoules_per_mole return EnergyTerms - +def evaluate_potential(Sim): + # Getting potential energy, taking RPMD into account(average over all copy of systems). If running classical simulation, it does nothing more than calling the getPotential energy function + if isinstance(Sim.integrator, RPMDIntegrator): + PE=0.0*kilojoule/mole + for i in range(Sim.integrator.getNumCopies()): + PE=PE+Sim.integrator.getState(i,getEnergy=True).getPotentialEnergy()/(Sim.integrator.getNumCopies()) + return PE + else: + return Sim.context.getState(getEnergy=True).getPotentialEnergy() +def evaluate_kinetic(Sim): + # It gets kinetic energy for classical simulation, or kinetic energy in classical sense(i.e. only dep. on how fast particles move) in RPMD simulation(sample to thermostat T*# of beads). It's NOT a quantum kinetic energy estimator. + if isinstance(Sim.integrator, RPMDIntegrator): + KE=0.0*kilojoule/mole + for i in range(Sim.integrator.getNumCopies()): + KE=KE+Sim.integrator.getState(i,getEnergy=True).getKineticEnergy()/(Sim.integrator.getNumCopies()) + return KE + else: + return Sim.context.getState(getEnergy=True).getKineticEnergy() +def primitive_kinetic(Sim): + # This is primitive quantum kinetic energy estimator for RPMD simulation. Return classical KE in classical simulation. + if isinstance(Sim.integrator,RPMDintegrator): + priKE=0.0*kilojoule/mole + hbar=0.0635078*nanometer**2*dalton/picosecond + for i in range(Sim.integrator.getNumCopies()): + priKE=priKE+Sim.integrator.getState(i,getEnergy=True).getKineticEnergy()/(Sim.integrator.getNumCopies()) #First term in primitive KE + mass_matrix=[] + for i in range(Sim.system.getNumParticles()): + mass_matrix.append(Sim.system.getParticleMass(i)) + mass_matrix=np.array(mass_matrix) + for i in range(Sim.integrator.getNumCopies()): + j=(i+1)%(Sim.integrator.getNumCopies()) + beaddif=np.array(Sim.integrator.getState(j,getPositions=True).getPositions())-np.array(Sim.integrator.getState(i,getPositions=True).getPositions())#calculate difference between ith and i+1th bead + priKE=priKE-np.sum(((beaddif*beaddif).sum(axis=1))*mass_matrix*(kB**2*Sim.integrator.getTemperature()**2*Sim.integrator.getNumCopies()/(2.0*hbar**2)))#2nd term in primitive estimator + return priKE + else: + return Sim.context.getState(getEnergy=True).getKineticEnergy() +def centroid_kinetic(Sim): + # This is centroid quantum kinetic energy estimator for RPMD simulation. Return classical KE in classical simulation. + if isinstance(Sim.integrator,RPMDintegrator): + cenKE=0.0*kilojoule/mole + for i in range(Sim.integrator.getNumCopies()): + cenKE=cenKE+Sim.integrator.getState(i,getEnergy=True).getKineticEnergy()/(Sim.integrator.getNumCopies()**2)#First term in centroid KE + centroid=np.array([[0.0*nanometer,0.0*nanometer,0.0*nanometer]]*Sim.system.getNumParticles()) + for i in range(Sim.integrator.getNumCopies()): + centroid=centroid+np.array(Sim.integrator.getState(i,getPositions=True).getPositions())/Sim.integrator.getNumCopies()#Calculate centroid of ring polymer + for i in range(Sim.integrator.getNumCopies()):#2nd term of centroid KE + dif=np.array(Sim.integrator.getState(i,getPositions=True).getPositions())-centroid + der=-1.0*np.array(Sim.integrator.getState(i,getForces=True).getForces()) + cenKE=cenKE+np.sum((dif*der).sum(axis=1))*0.5/Sim.integrator.getNumCopies() + return cenKE + else: + return Sim.context.getState(getEnergy=True).getKineticEnergy() def get_multipoles(simulation,q=None,mass=None,positions=None,rmcom=True): """Return the current multipole moments in Debye and Buckingham units. """ dx = 0.0 From 83817f3009d5eb60c91712e58780f5a88b23b8ad Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Tue, 8 Apr 2014 17:45:14 -0700 Subject: [PATCH 20/34] Replace code that calculate quantum kinetic energy by function --- src/openmmio.py | 80 +++++++------------------------------------------ 1 file changed, 10 insertions(+), 70 deletions(-) diff --git a/src/openmmio.py b/src/openmmio.py index 228fa6e44..de0db7063 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -1114,41 +1114,11 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, self.simulation.step(nsave) state = self.simulation.context.getState(getEnergy=True,getPositions=True,getVelocities=False,getForces=False) -##### Virial Energy estimator - if self.tdiv>=2: - hbar=0.0635078*nanometer**2*dalton/picosecond - pimdstate=[] - centroid=np.array([[0.0*nanometer,0.0*nanometer,0.0*nanometer]]*self.system.getNumParticles()) - kinetic=0.0*kilojoule/mole - cent_kinetic=0.0*kilojoule/mole - pri_kinetic=0.0*kilojoule/mole - potential=0.0*kilojoule/mole - integrator=self.simulation.context.getIntegrator() - for i in range(self.tdiv): - pimdstate.append(integrator.getState(i,getEnergy=True,getPositions=True,getForces=True,groups=-1)) - for i in range(self.tdiv): - centroid=centroid+np.array(pimdstate[i].getPositions())/self.tdiv - potential=potential+pimdstate[i].getPotentialEnergy()/self.tdiv#potential energy - cent_kinetic=cent_kinetic+pimdstate[i].getKineticEnergy()/(self.tdiv*self.tdiv)#1st term in centroid energy estimator - pri_kinetic=pri_kinetic+pimdstate[i].getKineticEnergy()/self.tdiv#1st term in primitive kinetic energy estimator - kinetic=kinetic+pimdstate[i].getKineticEnergy()/(self.tdiv*self.tdiv)#the kinetic energy that average to thermostat temperature - for i in range(self.tdiv):#2nd term in centroid kinetic energy estimator - dif=np.array(pimdstate[i].getPositions())-centroid - der=-1.0*np.array(pimdstate[i].getForces()) - cent_kinetic=cent_kinetic+np.sum((dif*der).sum(axis=1))*0.5/self.tdiv - mass_matrix=[] - for i in range(self.system.getNumParticles()): - mass_matrix.append(self.system.getParticleMass(i)) - mass_matrix=np.array(mass_matrix) - for i in range(self.tdiv):#2nd term in primitive kinetic energy estimator - j=(i+1)%(self.tdiv) - beaddif=np.array(pimdstate[j].getPositions())-np.array(pimdstate[i].getPositions())#difference in position between ith and i+1th bead - pri_kinetic=pri_kinetic-np.sum(((beaddif*beaddif).sum(axis=1))*mass_matrix*(kB**2*integrator.getTemperature()**2*self.tdiv)/(2.0*hbar**2)) - else: - kinetic = state.getKineticEnergy()/self.tdiv - potential = state.getPotentialEnergy() - cent_kinetic=0.0*kilojoule/mole - pri_kinetic=0.0*kilojoule/mole#Don't use centroid or primitive kinetic energy estimator for non-RPMD simulation +##### Energy Data + kinetic=evaluate_potential(self.simulation) + potential=evaluate_kinetic(self.simulation) + pri_kinetic=primitive_kinetic(self.simulation) + cen_kinetic=centroid_kinetic(self.simulation) ##### if self.pbc: box_vectors = state.getPeriodicBoxVectors() @@ -1179,41 +1149,11 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, if iteration >= 0: self.simulation.step(nsave) # Compute properties. state = self.simulation.context.getState(getEnergy=True,getPositions=True,getVelocities=False,getForces=False) -#####Virial Energy estimator - if self.tdiv>=2: - hbar=0.0635078*nanometer**2*dalton/picosecond - pimdstate=[] - centroid=np.array([[0.0*nanometer,0.0*nanometer,0.0*nanometer]]*self.system.getNumParticles()) - kinetic=0.0*kilojoule/mole - cent_kinetic=0.0*kilojoule/mole - pri_kinetic=0.0*kilojoule/mole - potential=0.0*kilojoule/mole - integrator=self.simulation.context.getIntegrator() - for i in range(self.tdiv): - pimdstate.append(integrator.getState(i,getEnergy=True,getPositions=True,getForces=True,groups=-1)) - for i in range(self.tdiv): - centroid=centroid+np.array(pimdstate[i].getPositions())/self.tdiv - potential=potential+pimdstate[i].getPotentialEnergy()/self.tdiv#potential energy - cent_kinetic=cent_kinetic+pimdstate[i].getKineticEnergy()/(self.tdiv*self.tdiv)#1st term in centroid energy estimator - pri_kinetic=pri_kinetic+pimdstate[i].getKineticEnergy()/self.tdiv#1st term in primitive kinetic energy estimator - kinetic=kinetic+pimdstate[i].getKineticEnergy()/(self.tdiv*self.tdiv)#the kinetic energy that average to thermostat temperature - for i in range(self.tdiv):#2nd term in centroid kinetic energy estimator - dif=np.array(pimdstate[i].getPositions())-centroid - der=-1.0*np.array(pimdstate[i].getForces()) - cent_kinetic=cent_kinetic+np.sum((dif*der).sum(axis=1))*0.5/self.tdiv - mass_matrix=[] - for i in range(self.system.getNumParticles()): - mass_matrix.append(self.system.getParticleMass(i)) - mass_matrix=np.array(mass_matrix) - for i in range(self.tdiv):#2nd term in primitive kinetic energy estimator - j=(i+1)%(self.tdiv) - beaddif=np.array(pimdstate[j].getPositions())-np.array(pimdstate[i].getPositions())#difference in position between ith and i+1th bead - pri_kinetic=pri_kinetic-np.sum(((beaddif*beaddif).sum(axis=1))*mass_matrix*(kB**2*integrator.getTemperature()**2*self.tdiv)/(2.0*hbar**2)) - else: - kinetic = state.getKineticEnergy()/self.tdiv - potential = state.getPotentialEnergy() - cent_kinetic=0.0*kilojoule/mole - pri_kinetic=0.0*kilojoule/mole#Don't use centroid or primitive kinetic energy estimator for non-RPMD simulation +##### Energy Data + kinetic=evaluate_potential(self.simulation) + potential=evaluate_kinetic(self.simulation) + pri_kinetic=primitive_kinetic(self.simulation) + cen_kinetic=centroid_kinetic(self.simulation) ##### kinetic_temperature = 2.0 * kinetic / kB / self.ndof if self.pbc: From fc50a07bad59cf3d369481e6b15aef6106e2eea2 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Wed, 9 Apr 2014 11:33:04 -0700 Subject: [PATCH 21/34] Fix a few RPMD related bug --- src/openmmio.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/openmmio.py b/src/openmmio.py index de0db7063..c1a44567c 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -41,6 +41,8 @@ def energy_components(Sim, verbose=False): for i in range(Sim.system.getNumForces()): EnergyTerms[Sim.system.getForce(i).__class__.__name__] = Sim.context.getState(getEnergy=True,groups=2**i).getPotentialEnergy() / kilojoules_per_mole return EnergyTerms +#def centroid_position(Sim): + # def evaluate_potential(Sim): # Getting potential energy, taking RPMD into account(average over all copy of systems). If running classical simulation, it does nothing more than calling the getPotential energy function if isinstance(Sim.integrator, RPMDIntegrator): @@ -61,9 +63,10 @@ def evaluate_kinetic(Sim): return Sim.context.getState(getEnergy=True).getKineticEnergy() def primitive_kinetic(Sim): # This is primitive quantum kinetic energy estimator for RPMD simulation. Return classical KE in classical simulation. - if isinstance(Sim.integrator,RPMDintegrator): + if isinstance(Sim.integrator,RPMDIntegrator): priKE=0.0*kilojoule/mole hbar=0.0635078*nanometer**2*dalton/picosecond + kb=0.00831446*nanometer**2*dalton/(picosecond**2*kelvin) for i in range(Sim.integrator.getNumCopies()): priKE=priKE+Sim.integrator.getState(i,getEnergy=True).getKineticEnergy()/(Sim.integrator.getNumCopies()) #First term in primitive KE mass_matrix=[] @@ -73,13 +76,13 @@ def primitive_kinetic(Sim): for i in range(Sim.integrator.getNumCopies()): j=(i+1)%(Sim.integrator.getNumCopies()) beaddif=np.array(Sim.integrator.getState(j,getPositions=True).getPositions())-np.array(Sim.integrator.getState(i,getPositions=True).getPositions())#calculate difference between ith and i+1th bead - priKE=priKE-np.sum(((beaddif*beaddif).sum(axis=1))*mass_matrix*(kB**2*Sim.integrator.getTemperature()**2*Sim.integrator.getNumCopies()/(2.0*hbar**2)))#2nd term in primitive estimator + priKE=priKE-np.sum(((beaddif*beaddif).sum(axis=1))*mass_matrix*(kb**2*Sim.integrator.getTemperature()**2*Sim.integrator.getNumCopies()/(2.0*hbar**2)))#2nd term in primitive estimator return priKE else: return Sim.context.getState(getEnergy=True).getKineticEnergy() def centroid_kinetic(Sim): # This is centroid quantum kinetic energy estimator for RPMD simulation. Return classical KE in classical simulation. - if isinstance(Sim.integrator,RPMDintegrator): + if isinstance(Sim.integrator,RPMDIntegrator): cenKE=0.0*kilojoule/mole for i in range(Sim.integrator.getNumCopies()): cenKE=cenKE+Sim.integrator.getState(i,getEnergy=True).getKineticEnergy()/(Sim.integrator.getNumCopies()**2)#First term in centroid KE From f1a5cdebfd4573c3e900b9bed09070c5c8fea655 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Wed, 9 Apr 2014 11:51:44 -0700 Subject: [PATCH 22/34] Adding centroid position calculator --- src/openmmio.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/openmmio.py b/src/openmmio.py index c1a44567c..28de9caa5 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -41,8 +41,15 @@ def energy_components(Sim, verbose=False): for i in range(Sim.system.getNumForces()): EnergyTerms[Sim.system.getForce(i).__class__.__name__] = Sim.context.getState(getEnergy=True,groups=2**i).getPotentialEnergy() / kilojoules_per_mole return EnergyTerms -#def centroid_position(Sim): - # +def centroid_position(Sim): + #Calculate RPMD centroid position in RPMD simulation(May not be what you actually need!), no difference from simply call getPosition in Classical simulation + if isinstance(Sim.integrator, RPMDIntegrator): + centroid=np.array([[0.0*nanometer,0.0*nanometer,0.0*nanometer]]*Sim.system.getNumParticles()) + for i in range(Sim.integrator.getNumCopies()): + centroid=centroid+np.array(Sim.integrator.getState(i,getPositions=True).getPositions())/Sim.integrator.getNumCopies()#Calculate centroid of ring polymer + return centroid + else: + return Sim.context.getState(getPositions=True).getPositions() def evaluate_potential(Sim): # Getting potential energy, taking RPMD into account(average over all copy of systems). If running classical simulation, it does nothing more than calling the getPotential energy function if isinstance(Sim.integrator, RPMDIntegrator): From 9cdf137f8661d190df7b569ce3aa80a5d4469e53 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Wed, 9 Apr 2014 13:56:38 -0700 Subject: [PATCH 23/34] Change nothing, just test --- src/openmmio.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/openmmio.py b/src/openmmio.py index 28de9caa5..a751c37d6 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -42,7 +42,7 @@ def energy_components(Sim, verbose=False): EnergyTerms[Sim.system.getForce(i).__class__.__name__] = Sim.context.getState(getEnergy=True,groups=2**i).getPotentialEnergy() / kilojoules_per_mole return EnergyTerms def centroid_position(Sim): - #Calculate RPMD centroid position in RPMD simulation(May not be what you actually need!), no difference from simply call getPosition in Classical simulation + ##Calculate RPMD centroid position in RPMD simulation(May not be what you actually need!), no difference from simply call getPosition in Classical simulation if isinstance(Sim.integrator, RPMDIntegrator): centroid=np.array([[0.0*nanometer,0.0*nanometer,0.0*nanometer]]*Sim.system.getNumParticles()) for i in range(Sim.integrator.getNumCopies()): From 8b69702e92bd7d93e27084402d28334ae6a56720 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Wed, 9 Apr 2014 16:56:03 -0700 Subject: [PATCH 24/34] Fix some minor(but important!) bug. --- src/openmmio.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/openmmio.py b/src/openmmio.py index a751c37d6..7b672476c 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -159,6 +159,9 @@ def get_multipoles(simulation,q=None,mass=None,positions=None,rmcom=True): qzz -= tr/3 # This ordering has to do with the way TINKER prints it out. return [dx,dy,dz,qxx,qxy,qyy,qxz,qyz,qzz] +def get_multipoles_ave(simulation,q=None,mass=None,positions=None,rmcom=True) +# If RPMD simulation, average over all copies of systems, else(in classical simulation) just do the same thing get_multipoles do + def get_dipole(simulation,q=None,mass=None,positions=None): """Return the current dipole moment in Debye. @@ -1125,8 +1128,8 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, state = self.simulation.context.getState(getEnergy=True,getPositions=True,getVelocities=False,getForces=False) ##### Energy Data - kinetic=evaluate_potential(self.simulation) - potential=evaluate_kinetic(self.simulation) + kinetic=evaluate_kinetic(self.simulation) + potential=evaluate_potential(self.simulation) pri_kinetic=primitive_kinetic(self.simulation) cen_kinetic=centroid_kinetic(self.simulation) ##### @@ -1160,8 +1163,8 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, # Compute properties. state = self.simulation.context.getState(getEnergy=True,getPositions=True,getVelocities=False,getForces=False) ##### Energy Data - kinetic=evaluate_potential(self.simulation) - potential=evaluate_kinetic(self.simulation) + kinetic=evaluate_kinetic(self.simulation) + potential=evaluate_potential(self.simulation) pri_kinetic=primitive_kinetic(self.simulation) cen_kinetic=centroid_kinetic(self.simulation) ##### From 520ede4215061e245f633547691e196e1d333bb9 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Wed, 9 Apr 2014 16:58:59 -0700 Subject: [PATCH 25/34] Try again. --- src/openmmio.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/openmmio.py b/src/openmmio.py index 7b672476c..0797bd85d 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -159,7 +159,7 @@ def get_multipoles(simulation,q=None,mass=None,positions=None,rmcom=True): qzz -= tr/3 # This ordering has to do with the way TINKER prints it out. return [dx,dy,dz,qxx,qxy,qyy,qxz,qyz,qzz] -def get_multipoles_ave(simulation,q=None,mass=None,positions=None,rmcom=True) +#def get_multipoles_ave(simulation,q=None,mass=None,positions=None,rmcom=True) # If RPMD simulation, average over all copies of systems, else(in classical simulation) just do the same thing get_multipoles do From 8a7d23db70101fcaf007228ee866d92c0b091f61 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Thu, 10 Apr 2014 11:53:02 -0700 Subject: [PATCH 26/34] RPMD Virtual Site --- src/openmmio.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/openmmio.py b/src/openmmio.py index 0797bd85d..642822926 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -845,11 +845,15 @@ def set_positions(self, shot=0, traj=None): self.simulation.context.setPositions(self.xyz_omms[shot][0]) self.simulation.context.computeVirtualSites() if self.tdiv>1: + integrator=self.simulation.context.getIntegrator() for i in range(self.tdiv): - integrator=self.simulation.context.getIntegrator() integrator.setPositions(i,self.xyz_omms[shot][0]) - ##compute virtual site method?^-^ - + for i in range(self.tdiv) + temp_position=integrator.getState(i,getPositions=true).getPositions() + self.simulation.context.setPositions(temp_position) + self.simulation.context.computeVirtualSites() + position_with_virtual_site=self.simulation.context.getPositions() + integrator.setPositions(i,position_with_virtual_site) def compute_volume(self, box_vectors): """ Compute the total volume of an OpenMM system. """ [a,b,c] = box_vectors From 96bea441f16f49c7cd1af7c59abaae00a207e3e3 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Thu, 10 Apr 2014 11:55:58 -0700 Subject: [PATCH 27/34] change nothing --- src/openmmio.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/openmmio.py b/src/openmmio.py index 642822926..477c966c0 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -844,7 +844,7 @@ def set_positions(self, shot=0, traj=None): # self.simulation.context.setPositions(ResetVirtualSites_fast(self.xyz_omms[shot][0], self.vsinfo)) self.simulation.context.setPositions(self.xyz_omms[shot][0]) self.simulation.context.computeVirtualSites() - if self.tdiv>1: + if self.tdiv>1:#RPMDintegrator integrator=self.simulation.context.getIntegrator() for i in range(self.tdiv): integrator.setPositions(i,self.xyz_omms[shot][0]) From 67b6c75f661caea6e18967a5beef298f8ecf70c7 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Thu, 10 Apr 2014 11:57:45 -0700 Subject: [PATCH 28/34] fix syntax error --- src/openmmio.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/openmmio.py b/src/openmmio.py index 477c966c0..ccd467968 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -848,7 +848,7 @@ def set_positions(self, shot=0, traj=None): integrator=self.simulation.context.getIntegrator() for i in range(self.tdiv): integrator.setPositions(i,self.xyz_omms[shot][0]) - for i in range(self.tdiv) + for i in range(self.tdiv): temp_position=integrator.getState(i,getPositions=true).getPositions() self.simulation.context.setPositions(temp_position) self.simulation.context.computeVirtualSites() From b826b1df89066fce34bb52cc1d60bb41e42c7f0f Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Thu, 10 Apr 2014 15:47:17 -0700 Subject: [PATCH 29/34] Use centroid KE to calculate KE in RPMD --- src/openmmio.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/openmmio.py b/src/openmmio.py index ccd467968..9aaefe895 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -1198,7 +1198,10 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, Temps.append(kinetic_temperature / kelvin) Rhos.append(density.value_in_unit(kilogram / meter**3)) Potentials.append(potential / kilojoules_per_mole) - Kinetics.append(kinetic / kilojoules_per_mole) + if self.tdiv>1: + Kinetics.append(cen_kinetic / kilojoules_per_mole)#If RPMD, we use centroid estimator to calculate quantum kinetic energy + else: + Kinetics.append(kinetic / kilojoules_per_mole) Volumes.append(volume / nanometer**3) Dips.append(get_dipole(self.simulation,positions=self.xyz_omms[-1][0])) Rhos = np.array(Rhos) From cc9abba7fa61f162d49858a8f8c132366453af1b Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Mon, 5 May 2014 14:51:23 -0700 Subject: [PATCH 30/34] Edit the fractionation ratio energy function --- src/data/npt.py | 6 +++- src/openmmio.py | 85 +++++++++++++++++++++++++++++++++---------------- 2 files changed, 63 insertions(+), 28 deletions(-) diff --git a/src/data/npt.py b/src/data/npt.py index 2035ec4ae..62a67b983 100755 --- a/src/data/npt.py +++ b/src/data/npt.py @@ -420,7 +420,7 @@ def main(): Volumes = prop_return['Volumes'] Dips = prop_return['Dips'] EDA = prop_return['Ecomps'] - + getone = prop_return['One'] # Create a bunch of physical constants. # Energies are in kJ/mol # Lengths are in nanometers. @@ -492,6 +492,10 @@ def main(): # Condensed phase properties and derivatives. # #==============================================# + #---- + # Just get one + #---- + #---- # Density #---- diff --git a/src/openmmio.py b/src/openmmio.py index 9aaefe895..8592d2b8c 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -34,6 +34,26 @@ except: pass +def H_spring_energy(Sim): # Spring energy in RPMD simulation,only count H atom, it can be used to calculate the HD fractionation ratio. If classical simulation, return 0. + if isinstance(Sim.integrator,RPMDIntegrator): + springE=0.0*kilojoule/mole + hbar=0.0635078*nanometer**2*dalton/picosecond + kb=0.00831446*nanometer**2*dalton/(picosecond**2*kelvin) + mass_matrix=[] + for i in range(Sim.system.getNumParticles()): + if (Sim.system.getParticleMass(i)<1.02*dalton):#Only count Hydrogen + mass_matrix.append(Sim.system.getParticleMass(i)) + else + mass_matrix.append(0.0*dalton) + mass_matrix=np.array(mass_matrix) + for i in range(Sim.integrator.getNumCopies()): + j=(i+1)%(Sim.integrator.getNumCopies()) + beaddif=np.array(Sim.integrator.getState(j,getPositions=True).getPositions())-np.array(Sim.integrator.getState(i,getPositions=True).getPositions())#calculate difference between ith and i+1th bead + springE=springE+np.sum(((beaddif*beaddif).sum(axis=1))*mass_matrix*(kb**2*Sim.integrator.getTemperature()**2*Sim.integrator.getNumCopies()/(2.0*hbar**2)))#spring energy + return springE + else: + return Sim.context.getState(getEnergy=True).getKineticEnergy() + def energy_components(Sim, verbose=False): # Before using EnergyComponents, make sure each Force is set to a different group. EnergyTerms = OrderedDict() @@ -51,11 +71,13 @@ def centroid_position(Sim): else: return Sim.context.getState(getPositions=True).getPositions() def evaluate_potential(Sim): - # Getting potential energy, taking RPMD into account(average over all copy of systems). If running classical simulation, it does nothing more than calling the getPotential energy function + # Getting potential energy, taking RPMD into account(average over 4 copies of systems). If running classical simulation, it does nothing more than calling the getPotential energy function if isinstance(Sim.integrator, RPMDIntegrator): PE=0.0*kilojoule/mole - for i in range(Sim.integrator.getNumCopies()): - PE=PE+Sim.integrator.getState(i,getEnergy=True).getPotentialEnergy()/(Sim.integrator.getNumCopies()) + nob=Sim.integrator.getNumCopies() + hello=int(nob/4) + for i in range(0,nob,hello): + PE=PE+Sim.integrator.getState(i,getEnergy=True).getPotentialEnergy()/4.0 return PE else: return Sim.context.getState(getEnergy=True).getPotentialEnergy() @@ -63,8 +85,10 @@ def evaluate_kinetic(Sim): # It gets kinetic energy for classical simulation, or kinetic energy in classical sense(i.e. only dep. on how fast particles move) in RPMD simulation(sample to thermostat T*# of beads). It's NOT a quantum kinetic energy estimator. if isinstance(Sim.integrator, RPMDIntegrator): KE=0.0*kilojoule/mole - for i in range(Sim.integrator.getNumCopies()): - KE=KE+Sim.integrator.getState(i,getEnergy=True).getKineticEnergy()/(Sim.integrator.getNumCopies()) + nob=Sim.integrator.getNumCopies() + hello=int(nob/4) + for i in range(0,nob,hello): + KE=KE+Sim.integrator.getState(i,getEnergy=True).getKineticEnergy()/4.0 return KE else: return Sim.context.getState(getEnergy=True).getKineticEnergy() @@ -91,15 +115,17 @@ def centroid_kinetic(Sim): # This is centroid quantum kinetic energy estimator for RPMD simulation. Return classical KE in classical simulation. if isinstance(Sim.integrator,RPMDIntegrator): cenKE=0.0*kilojoule/mole - for i in range(Sim.integrator.getNumCopies()): - cenKE=cenKE+Sim.integrator.getState(i,getEnergy=True).getKineticEnergy()/(Sim.integrator.getNumCopies()**2)#First term in centroid KE + nob=Sim.integrator.getNumCopies() + hello=int(nob/4) + for i in range(0,nob,hello): + cenKE=cenKE+Sim.integrator.getState(i,getEnergy=True).getKineticEnergy()/(nob*4.0)#First term in centroid KE centroid=np.array([[0.0*nanometer,0.0*nanometer,0.0*nanometer]]*Sim.system.getNumParticles()) - for i in range(Sim.integrator.getNumCopies()): - centroid=centroid+np.array(Sim.integrator.getState(i,getPositions=True).getPositions())/Sim.integrator.getNumCopies()#Calculate centroid of ring polymer - for i in range(Sim.integrator.getNumCopies()):#2nd term of centroid KE + for i in range(0,nob,hello): + centroid=centroid+np.array(Sim.integrator.getState(i,getPositions=True).getPositions())/4.0#Calculate centroid of ring polymer + for i in range(0,nob,hello):#2nd term of centroid KE dif=np.array(Sim.integrator.getState(i,getPositions=True).getPositions())-centroid der=-1.0*np.array(Sim.integrator.getState(i,getForces=True).getForces()) - cenKE=cenKE+np.sum((dif*der).sum(axis=1))*0.5/Sim.integrator.getNumCopies() + cenKE=cenKE+np.sum((dif*der).sum(axis=1))*0.5/4.0 return cenKE else: return Sim.context.getState(getEnergy=True).getKineticEnergy() @@ -780,6 +806,7 @@ def create_simulation(self, timestep=1.0, faststep=0.25, temperature=None, press ## Finally create the simulation object. self.simulation = Simulation(self.mod.topology, self.system, integrator, self.platform) + #fffff.write(XmlSerializer.serialize(self.simulation.system)) ## Print platform properties. # logger.info("I'm using the platform %s\n" % self.simulation.context.getPlatform().getName()) # printcool_dictionary({i:self.simulation.context.getPlatform().getPropertyValue(self.simulation.context,i) \ @@ -848,12 +875,12 @@ def set_positions(self, shot=0, traj=None): integrator=self.simulation.context.getIntegrator() for i in range(self.tdiv): integrator.setPositions(i,self.xyz_omms[shot][0]) - for i in range(self.tdiv): - temp_position=integrator.getState(i,getPositions=true).getPositions() - self.simulation.context.setPositions(temp_position) - self.simulation.context.computeVirtualSites() - position_with_virtual_site=self.simulation.context.getPositions() - integrator.setPositions(i,position_with_virtual_site) +# for i in range(self.tdiv): +# temp_position=integrator.getState(i,getPositions=True).getPositions() +# self.simulation.context.setPositions(temp_position) +# self.simulation.context.computeVirtualSites() +# position_with_virtual_site=self.simulation.context.getState(getPositions=True).getPositions() +# integrator.setPositions(i,position_with_virtual_site) def compute_volume(self, box_vectors): """ Compute the total volume of an OpenMM system. """ [a,b,c] = box_vectors @@ -1118,7 +1145,7 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, #========================# # Initialize velocities. self.simulation.context.setVelocitiesToTemperature(temperature*kelvin) - + #print(XmlSerializer.serialize(self.simulation.system)) # Equilibrate. if iequil > 0: if verbose: logger.info("Equilibrating...\n") @@ -1126,15 +1153,17 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, if verbose: logger.info("%6s %9s %9s %13s %10s %13s\n" % ("Iter.", "Time(ps)", "Temp(K)", "Epot(kJ/mol)", "Vol(nm^3)", "Rho(kg/m^3)")) else: if verbose: logger.info("%6s %9s %9s %13s\n" % ("Iter.", "Time(ps)", "Temp(K)", "Epot(kJ/mol)")) - for iteration in range(-1, iequil): + for iteration in range(-1 if self.tdiv == 1 else 0, iequil): if iteration >= 0: self.simulation.step(nsave) state = self.simulation.context.getState(getEnergy=True,getPositions=True,getVelocities=False,getForces=False) ##### Energy Data - kinetic=evaluate_kinetic(self.simulation) +# kinetic=state.getKineticEnergy()/self.tdiv +# potential=state.getPotentialEnergy() + kinetic=evaluate_kinetic(self.simulation)/self.tdiv potential=evaluate_potential(self.simulation) - pri_kinetic=primitive_kinetic(self.simulation) + #pri_kinetic=primitive_kinetic(self.simulation) cen_kinetic=centroid_kinetic(self.simulation) ##### if self.pbc: @@ -1161,15 +1190,17 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, if save_traj: self.simulation.reporters.append(PDBReporter('%s-md.pdb' % self.name, nsteps)) self.simulation.reporters.append(DCDReporter('%s-md.dcd' % self.name, nsave)) - for iteration in range(-1, isteps): + for iteration in range(-1 if self.tdiv == 1 else 0, isteps): # Propagate dynamics. if iteration >= 0: self.simulation.step(nsave) # Compute properties. state = self.simulation.context.getState(getEnergy=True,getPositions=True,getVelocities=False,getForces=False) ##### Energy Data - kinetic=evaluate_kinetic(self.simulation) +# kinetic=state.getKineticEnergy()/self.tdiv +# potential=state.getPotentialEnergy() + kinetic=evaluate_kinetic(self.simulation)/self.tdiv potential=evaluate_potential(self.simulation) - pri_kinetic=primitive_kinetic(self.simulation) + #pri_kinetic=primitive_kinetic(self.simulation) cen_kinetic=centroid_kinetic(self.simulation) ##### kinetic_temperature = 2.0 * kinetic / kB / self.ndof @@ -1199,9 +1230,9 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, Rhos.append(density.value_in_unit(kilogram / meter**3)) Potentials.append(potential / kilojoules_per_mole) if self.tdiv>1: - Kinetics.append(cen_kinetic / kilojoules_per_mole)#If RPMD, we use centroid estimator to calculate quantum kinetic energy + Kinetics.append(cen_kinetic / kilojoules_per_mole)#If RPMD, we use centroid estimator to calculate quantum kinetic energy else: - Kinetics.append(kinetic / kilojoules_per_mole) + Kinetics.append(kinetic / kilojoules_per_mole)#else just classical kinetic energy Volumes.append(volume / nanometer**3) Dips.append(get_dipole(self.simulation,positions=self.xyz_omms[-1][0])) Rhos = np.array(Rhos) @@ -1216,7 +1247,7 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, Ecomps["Total Energy"] = np.array(Potentials) + np.array(Kinetics) # Initialized property dictionary. prop_return = OrderedDict() - prop_return.update({'Rhos': Rhos, 'Potentials': Potentials, 'Kinetics': Kinetics, 'Volumes': Volumes, 'Dips': Dips, 'Ecomps': Ecomps}) + prop_return.update({'Rhos': Rhos, 'Potentials': Potentials, 'Kinetics': Kinetics, 'Volumes': Volumes, 'Dips': Dips, 'Ecomps': Ecomps, 'One': 1.0}) return prop_return class Liquid_OpenMM(Liquid): From 9a8950262be029e583b020981b4d10296ceb704e Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Mon, 5 May 2014 14:55:03 -0700 Subject: [PATCH 31/34] edit HD --- src/openmmio.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/openmmio.py b/src/openmmio.py index 8592d2b8c..1252b57bd 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -52,7 +52,7 @@ def H_spring_energy(Sim): # Spring energy in RPMD simulation,only count H atom, springE=springE+np.sum(((beaddif*beaddif).sum(axis=1))*mass_matrix*(kb**2*Sim.integrator.getTemperature()**2*Sim.integrator.getNumCopies()/(2.0*hbar**2)))#spring energy return springE else: - return Sim.context.getState(getEnergy=True).getKineticEnergy() + return 0.0*kilojoule/mole def energy_components(Sim, verbose=False): # Before using EnergyComponents, make sure each Force is set to a different group. From 5a427d9f74497ad40b00c4cb26f144f77d0740b0 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Mon, 5 May 2014 14:58:40 -0700 Subject: [PATCH 32/34] ":" --- src/openmmio.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/openmmio.py b/src/openmmio.py index 1252b57bd..a2b29e7d5 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -43,7 +43,7 @@ def H_spring_energy(Sim): # Spring energy in RPMD simulation,only count H atom, for i in range(Sim.system.getNumParticles()): if (Sim.system.getParticleMass(i)<1.02*dalton):#Only count Hydrogen mass_matrix.append(Sim.system.getParticleMass(i)) - else + else: mass_matrix.append(0.0*dalton) mass_matrix=np.array(mass_matrix) for i in range(Sim.integrator.getNumCopies()): From 09944082dd4cfb47dd0dae70043f6c1a04ddca35 Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Mon, 5 May 2014 15:13:31 -0700 Subject: [PATCH 33/34] undo last step --- src/openmmio.py | 38 +++++++++++++++++++------------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/src/openmmio.py b/src/openmmio.py index a2b29e7d5..346459f0b 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -34,25 +34,25 @@ except: pass -def H_spring_energy(Sim): # Spring energy in RPMD simulation,only count H atom, it can be used to calculate the HD fractionation ratio. If classical simulation, return 0. - if isinstance(Sim.integrator,RPMDIntegrator): - springE=0.0*kilojoule/mole - hbar=0.0635078*nanometer**2*dalton/picosecond - kb=0.00831446*nanometer**2*dalton/(picosecond**2*kelvin) - mass_matrix=[] - for i in range(Sim.system.getNumParticles()): - if (Sim.system.getParticleMass(i)<1.02*dalton):#Only count Hydrogen - mass_matrix.append(Sim.system.getParticleMass(i)) - else: - mass_matrix.append(0.0*dalton) - mass_matrix=np.array(mass_matrix) - for i in range(Sim.integrator.getNumCopies()): - j=(i+1)%(Sim.integrator.getNumCopies()) - beaddif=np.array(Sim.integrator.getState(j,getPositions=True).getPositions())-np.array(Sim.integrator.getState(i,getPositions=True).getPositions())#calculate difference between ith and i+1th bead - springE=springE+np.sum(((beaddif*beaddif).sum(axis=1))*mass_matrix*(kb**2*Sim.integrator.getTemperature()**2*Sim.integrator.getNumCopies()/(2.0*hbar**2)))#spring energy - return springE - else: - return 0.0*kilojoule/mole +#def H_spring_energy(Sim): # Spring energy in RPMD simulation,only count H atom, it can be used to calculate the HD fractionation ratio. If classical simulation, return 0. +# if isinstance(Sim.integrator,RPMDIntegrator): +# springE=0.0*kilojoule/mole +# hbar=0.0635078*nanometer**2*dalton/picosecond +# kb=0.00831446*nanometer**2*dalton/(picosecond**2*kelvin) +# mass_matrix=[] +# for i in range(Sim.system.getNumParticles()): +# if (Sim.system.getParticleMass(i)<1.02*dalton):#Only count Hydrogen +# mass_matrix.append(Sim.system.getParticleMass(i)) +# else: +# mass_matrix.append(0.0*dalton) +# mass_matrix=np.array(mass_matrix) +# for i in range(Sim.integrator.getNumCopies()): +# j=(i+1)%(Sim.integrator.getNumCopies()) +# beaddif=np.array(Sim.integrator.getState(j,getPositions=True).getPositions())-np.array(Sim.integrator.getState(i,getPositions=True).getPositions())#calculate difference between ith and i+1th bead +# springE=springE+np.sum(((beaddif*beaddif).sum(axis=1))*mass_matrix*(kb**2*Sim.integrator.getTemperature()**2*Sim.integrator.getNumCopies()/(2.0*hbar**2)))#spring energy +# return springE +# else: +# return 0.0*kilojoule/mole def energy_components(Sim, verbose=False): # Before using EnergyComponents, make sure each Force is set to a different group. From a80f321dacba2c8f14fc7a84531b40cca0a0e29e Mon Sep 17 00:00:00 2001 From: Frank Lin Date: Tue, 6 May 2014 13:09:08 -0700 Subject: [PATCH 34/34] undo --- src/data/npt.py | 2 +- src/openmmio.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/data/npt.py b/src/data/npt.py index 62a67b983..64d05b2ed 100755 --- a/src/data/npt.py +++ b/src/data/npt.py @@ -420,7 +420,7 @@ def main(): Volumes = prop_return['Volumes'] Dips = prop_return['Dips'] EDA = prop_return['Ecomps'] - getone = prop_return['One'] + #getone = prop_return['One'] # Create a bunch of physical constants. # Energies are in kJ/mol # Lengths are in nanometers. diff --git a/src/openmmio.py b/src/openmmio.py index 346459f0b..fdd708df6 100644 --- a/src/openmmio.py +++ b/src/openmmio.py @@ -1247,7 +1247,7 @@ def molecular_dynamics(self, nsteps, timestep, temperature=None, pressure=None, Ecomps["Total Energy"] = np.array(Potentials) + np.array(Kinetics) # Initialized property dictionary. prop_return = OrderedDict() - prop_return.update({'Rhos': Rhos, 'Potentials': Potentials, 'Kinetics': Kinetics, 'Volumes': Volumes, 'Dips': Dips, 'Ecomps': Ecomps, 'One': 1.0}) + prop_return.update({'Rhos': Rhos, 'Potentials': Potentials, 'Kinetics': Kinetics, 'Volumes': Volumes, 'Dips': Dips, 'Ecomps': Ecomps}) return prop_return class Liquid_OpenMM(Liquid):