search

protocaller / ProtoCaller

Full automation of relative protein-ligand binding free energy calculations in GROMACS
http://protocaller.readthedocs.io
GNU General Public License v3.0
43 stars 15 forks source link

LINCS warning and structure clash #27

Closed kexul closed 3 years ago

kexul commented 3 years ago

Hi, here is my transformation image

here is the .gro and .top file generate by ProtoCaller, morphs.zip

The minimization step seems fine, but, when running NVT, there are many LINCS warning, below are some of them: 

Step 9209, time 9.209 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000115, max 0.000551 (between atoms 47 and 58)
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length
     47     57   30.5    0.1097   0.1097      0.1097

Step 9210, time 9.21 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000074, max 0.000297 (between atoms 35 and 49)
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length
     47     57   31.5    0.1097   0.1097      0.1097

When I check the gro file generate after NVT Equilibration, atom 49 and atom 43 are too close, which I think should not happen in reality. image

Here is the .mdp file I used and .gro file generated after NVT Equilibration. files.zip

Any clue what happened? Thanks in advance!

kexul commented 3 years ago

I've consulted an expert, who told me that the .top file generated by protocaller had some missing dihedrals parameters. image

kexul commented 3 years ago

Besides, the dummy atom's name seems to be duplicated with previous atoms, would that be a problem ?

image

msuruzhon commented 3 years ago

It does seem that there are issues with the dummy atoms but it's not obvious what the nature of these is. I think the best course of action is to try the opposite perturbation which should result in unique atom names and try running these. It is also worth comparing the dihedrals from each separate ligand and make sure that they are correctly represented in the merged topology, but it is also possible that GAFF doesn't return any dihedral terms for this particular bond, so it is set to 0.

msuruzhon commented 3 years ago

I just looked at your mdp file. You are turning on the electrostatic interactions for ligand B before turning on the Lennard-Jones interactions, and this is most likely the reason for the behaviour you are observing.

kexul commented 3 years ago

I just looked at your mdp file. You are turning on the electrostatic interactions for ligand B before turning on the Lennard-Jones interactions, and this is most likely the reason for the behaviour you are observing.

Wow! I changed the order as you said, which did the trick. Thanks so much for your advice! By the way, the old order I followed was taken from the example of ProtoCaller, except increasing the number of lambda windows. Would you mind elaborate the difference between the cases? In which situation should we turning on the LJ before electrostatic?

tsenapathi commented 3 years ago

@kexul you are changing a small molecule to a larger one. Therefore you need to turn on the LJ before electrostatic; otherwise, charges will appear in the empty space which will lead to instabilities of the simulation.

kexul commented 3 years ago

@tsenapathi Got it, Thank you!