Energy decomposition of System using CHARMM36 FF
Posted: Mon Apr 11, 2016 2:33 am
I have been trying to do energy decomposition of a system running CHARMM36 FF on OpenMM.
I am using scripts provided by CHARMM-GUI to run the simulation because they have implemented CustomNonBondedForce to do force-based switching.
After creating the system and adding in the coordinates I ran the following script to get the energy decomposition:
The total potential energy is similar to the potential energy calculated by NAMD. But I am unable to reproduce the total potential energy from the energy decomposition. I get the following output by running the function (first line is the total potential energy):
1. Bond (Excludes Hbonds, using rigid HBonds)
2. Angle
3. Bond (Urey-Bradley term)
4. Dihedral
5. Improper
6. Cross-terms
7. Non-bonded (electrostatics)
8. Custom-nonbonded (van der Waals)
9. Center of mass motion remover
10. Custom-nonbonded (van der Waals)
11. Custom-bonded (van der Waals)
I have attached the python script that implements the force switching for the van der Waals.
Another thing that I noticed is that the number of dihedral angles reported by PeriodicTorsionForce is a lot more than the number report by NAMD/VMD and in my PSF file. Weirdly, when I run this system on AMBER (using CHAMBER to convert PSF) I get the exact same number as in OpenMM.
Thanks in advance!
I am using scripts provided by CHARMM-GUI to run the simulation because they have implemented CustomNonBondedForce to do force-based switching.
After creating the system and adding in the coordinates I ran the following script to get the energy decomposition:
Code: Select all
def getEnergyDecomposition(context, system):
print(context.getState(getEnergy=True).getPotentialEnergy())
forces = []
for i in range(system.getNumForces()):
current_force = system.getForce(i)
forces.append(current_force)
print(current_force, ':', context.getState(getEnergy=True, groups=2**i).getPotentialEnergy())
return forces
I have managed identify each force type to the FF terms:-1542535.43777 kJ/mol
<simtk.openmm.openmm.HarmonicBondForce; proxy of <Swig Object of type 'OpenMM::HarmonicBondForce *' at 0x7f5c26948060> > : 76593.1529571 kJ/mol
<simtk.openmm.openmm.HarmonicAngleForce; proxy of <Swig Object of type 'OpenMM::HarmonicAngleForce *' at 0x7f5c26948090> > : 137966.790553 kJ/mol
<simtk.openmm.openmm.HarmonicBondForce; proxy of <Swig Object of type 'OpenMM::HarmonicBondForce *' at 0x7f5c26948240> > : 150203.559384 kJ/mol
<simtk.openmm.openmm.PeriodicTorsionForce; proxy of <Swig Object of type 'OpenMM::PeriodicTorsionForce *' at 0x7f5c26948120> > : 45353.6148533 kJ/mol
<simtk.openmm.openmm.CustomTorsionForce; proxy of <Swig Object of type 'OpenMM::CustomTorsionForce *' at 0x7f5c269481b0> > : 8677.96254746 kJ/mol
<simtk.openmm.openmm.CMAPTorsionForce; proxy of <Swig Object of type 'OpenMM::CMAPTorsionForce *' at 0x7f5c26948180> > : 2961.66222989 kJ/mol
<simtk.openmm.openmm.NonbondedForce; proxy of <Swig Object of type 'OpenMM::NonbondedForce *' at 0x7f5c269481e0> > : -1935585.15758 kJ/mol
<simtk.openmm.openmm.CustomNonbondedForce; proxy of <Swig Object of type 'OpenMM::CustomNonbondedForce *' at 0x7f5c26948210> > : 4784.50361701 kJ/mol
<simtk.openmm.openmm.CMMotionRemover; proxy of <Swig Object of type 'OpenMM::CMMotionRemover *' at 0x7f5c26948270> > : 4784.50361701 kJ/mol
<simtk.openmm.openmm.CustomNonbondedForce; proxy of <Swig Object of type 'OpenMM::CustomNonbondedForce *' at 0x7f5c269482a0> > : 4784.50361701 kJ/mol
<simtk.openmm.openmm.CustomBondForce; proxy of <Swig Object of type 'OpenMM::CustomBondForce *' at 0x7f5c269482d0> > : 4784.50361701 kJ/mol
1. Bond (Excludes Hbonds, using rigid HBonds)
2. Angle
3. Bond (Urey-Bradley term)
4. Dihedral
5. Improper
6. Cross-terms
7. Non-bonded (electrostatics)
8. Custom-nonbonded (van der Waals)
9. Center of mass motion remover
10. Custom-nonbonded (van der Waals)
11. Custom-bonded (van der Waals)
I have attached the python script that implements the force switching for the van der Waals.
Another thing that I noticed is that the number of dihedral angles reported by PeriodicTorsionForce is a lot more than the number report by NAMD/VMD and in my PSF file. Weirdly, when I run this system on AMBER (using CHAMBER to convert PSF) I get the exact same number as in OpenMM.
Thanks in advance!