The NAMD scripts calculate the energies of a solvated protein-ligand QMMM system, including pair interaction energies. Software: NAMD 2.14 / MOPAC 2016. The QM part of the system is a ligand molecule (segid LIG). The MM part of the system is the protein+TIP3 water (with the total charge of 0, no counterions need to be included). Main files: Protein_Ligand_solvated.psf / Protein_Ligand_solvated.pdb : the protein-ligand system SansLigand.psf / SansLigand.pdb : same with the ligand removed Forcefield: CHARMM 36 + the ligand forcefield generated by CGenFF Test coordinates /velocities: production.coor / production.vel and production_sansligand.coor / production_sansligand.vel (final snapshots of a ~150ps MD run) Problem: The non-bonded, electrostatic and van der Waals energies dont' add up, especially the vdW energy that shoots up to crazy values during the pair interaction calculations. Tests: Test run E_elec E_vdw Comment E_test_no_PairInt.{namd,log} -932364.2484 -427.4665 The total system energy E_test_PairInt.{namd,log} -118280.8276 48861334.1618 QM part/MM part interaction energy E_test_PairInt_11.{namd.log} -118339.2930 48861334.1618 QM part/QM part interaction energy E_test_PairInt_00.{namd,log} -932294.1036 -427.4665 MM part/MM part interaction energy E_test_no_PairInt_noLigand.{namd,log} -932294.1034 -382.2416 MM part total energy with ligand removed Comment: The total system energy looks all right, in line with what was printed during the simulation. However, the QM-MM interaction energy makes no sense: why would the vdW component shoot up so much? Neither does the energy of the QM part interaction with itself: the crazy value shows again whereas the vdW energy should simply be zero in this case. NAMD doesn't calculate the vdW terms betwen the QM atoms, does it?? The pair interaction energy of the MM part with itself looks all right. However, why does the vdW energy have the same value as when I calculated the interaction of the whole system with itself?? Shouldn't it be E_vdw(MM-MM) = E_vdw(All-All) - E_vdw(QM-MM) ??? Furthermore, if I simply remove the ligand from the system and repeat the energy calculations, I am reasonable numbers again. However, the vdW energy is significantly different from that calculated for the MM-MM pair interaction. Why? shouldn't the two be the same? I further checked if it is the fixed atoms option that screws the calculations (I am fixing many atoms to speed up the QMMM run). Apparently, no: E_elec E_vdw E_test_no_PairInt_noFixed_noLigand.{namd,log} -293957.0182 23963.9304 E_test_PairInt_noFixed.{namd,log} -3827529.7215 124983177.2010 While the fixed atoms do contribute some positive vdW energy to the system, but it's not even comparable with the crazy vdW energy. The crazy energy does show up in pair interaction calculations again even if no atoms in the system are fixed.