ERROR: LINCS WARNINGS in WT-MetaD simulation

[Subarna Sasmal]

Hi everyone, 

I have been trying to run some WT-MetaD simulations using HLDA CV, for the vinculin system. The problem is, I am getting some LINCS WARNING  which I am not able to solve and that is why the simulations are crashing in between (sometimes just after a few ns).

sample errors -

Step 1781124, time 3562.25 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 9091.672852, max 62963.968750 (between atoms 951 and 952)
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length
    919    920   90.0    0.1080 461.0422      0.1080
    951    952   90.0    0.1080 6800.2168      0.1080
    973    974   90.0    0.1080 232.8021      0.1080
   1936   1937   90.0    0.1080 2735.6399      0.1080

Things I tried -

1. I have well equilibrated the system and  was able to run normal MD (200ns) using GROMACS with no interruption at all.
2.  I have tried different values of parameters for my WT-MetaD simulations.
3.  I checked if the bias getting added with time is too much at some point for some reason. But I didn't find anything like that.
4. Energy data from the simulations were monitored, nothing seemed abnormal.
5. Tried different versions of GROMACS-PLUMED.

But nothing could actually solve the issue.

I am attaching one sample PLUMED input file and one slurm output file for reference along with this email. Please let me know if you have any idea about how to deal with this issue. 

Thanks in advance.

Best,

Subarna Sasmal

NYU

["ROPÓN-PALACIOS G."]

Estimated users,

I'm new to plumed and I really don't have much idea how to start with any particular point.

My question is:

I have a protein system (prot-prot complex) and I want to measure mechanical properties, for this I want to perform SMD, but I don't know how to do this with plumed, and which version to use where it is implemented, as well as to perform a harmonic restriction of the coordinates of one of the proteins and finally how to perform umbrella sampling, with it.

Also if you can give me suggestions on how to correctly pass group of atoms.

If you can give me examples of the code I would appreciate it.

A thousand apologies, for the request for help, so trivial, for many.

Geo.

[Michele Invernizzi]

Hi Subarna Sasmal,

Does the error appear also if you remove the STRIDE=2 keyword from metad and the walls?

Best,

Michele

--
You received this message because you are subscribed to the Google Groups "PLUMED users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to plumed-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/eb10fbab-438b-4849-ae35-6f8f7ba25da1n%40googlegroups.com.

[Michele Invernizzi]

Hi Geo,

You should check out the PLUMED Masterclass tutorials https://www.plumed.org/doc-v2.7/user-doc/html/masterclass-21-1.html

There is one also on umbrella sampling.

And you can find the related theory lectures on YouTube https://youtube.com/channel/UCNboSvPixTVfsX91Fqli-8Q

Best,

Michele

--

You received this message because you are subscribed to the Google Groups "PLUMED users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to plumed-users...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/EDD2AA08-A0BD-48AA-98F8-2109010BA458%40gmail.com.

["ROPÓN-PALACIOS G."]

Thanks so much Michelle for advice. 

Best, 

Geo. 

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CAGBhpieHr3MsJazgU%3DY5Tmr0Eabc0Q_qiEsHiDxy4aXhzbgQDQ%40mail.gmail.com.

["ROPÓN-PALACIOS G."]

Dear users I’m using following code for performing SMD:

MOLINFO STRUCTURE=prot.pdb PYTHON_BIN=python3.7 

fix:    GROUP ATOMS={@mda: {segued PROA and name CA}}

Pull:  GROUP ATOMS={@mda: {segued PROB and name CA}}

d1: DISTANCE ATOMS=fix, pull COMPONENTS 

MOVINGRESTRAINT …

ARG=d1.z 

STEP=0    AT0=3.45 KAPPA0=1000.0 

STREP1=2500   AT1=4.0 KAPPA1=1000.0

…

PRINT ARG=d1.z FILE=D1.z STRIDE=100

I would like to have an orientation of if I am getting on with the syntax and resolve some doubts.

1. How do I assign or calculate the pulling speed?

2. As within this code I restraint the alpha carbons of the protein to be fixed (PROA)

3. In the output, as I print only the data of the distance in Z.

Greetings,

Geo.

PS: I'm just learning from this, a thousand apologies for the trivial questions

On Nov 1, 2021, at 03:26, Michele Invernizzi <michele.i...@fu-berlin.de> wrote:

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CAGBhpieHr3MsJazgU%3DY5Tmr0Eabc0Q_qiEsHiDxy4aXhzbgQDQ%40mail.gmail.com.

[Gareth Tribello]

Hello

You calculate the pulling speed by using distance over time so in your case this is:

(4 - 3.45) / (2500*delta)

where delta is your simulation timestep.

I am not sure how to answer your other two doubts as what you have written are statements and not questions.

I hope this helps

Gareth

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/AA63419A-429E-466F-9085-0EC82B460F28%40gmail.com.

["ROPÓN-PALACIOS G."]

Dear Tribello, 

I like : 

1. Restraint CHAIN A of my protein, why it is fixed during pulling , how make it? 

2.  What is the right definition of AT values into SMD bias? 

I have change my code to it below, using a virtual atoms as reference to calculate distance of pulling , in it point , 

I want know if STEP1 value is same that it steps to run into MDP config file of gromacs to run!

3. AT and KAPPA values is in nanometers and kj/mol*nm^2 respectively? 

4. What is unit of velocity of pulling, nm/ps? 

MOLINFO STRUCTURE=prot.pdb PYTHON_BIN=python3.7 

fix:    FIXEDATOM ATOMS=0,0,0

Pull:  CENTER ATOMS={@mda: {segued PROB and name CA}}

d1: DISTANCE ATOMS=fix, pull COMPONENTS 

MOVINGRESTRAINT …

ARG=d1.z 

STEP=0     AT0=3.45 KAPPA0=1000.0 

STEP1=2500   AT1=4.0 KAPPA1=1000.0

…

PRINT ARG=d1.z FILE=D1.z STRIDE=100

Best, 

Geo. 

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/6ED07442-5D8A-4ADB-8866-40D5C9549A10%40gmail.com.

[Subarna Sasmal]

Hi Michele,

Thanks for your reply.

I have tried running the simulations by removing the STRIDE=2 from the metad and walls. But it didn't improve anything because simulations are still crashing. 

I still got LINCS WARNINGS like before and also a few more errors which say the system is looking for values outside the grid range for HLDA CV. I checked the COLVAR files to see if the system is actually crossing the range during the simulations but found that it didn't reach close to the boundaries, so it didn't make sense.

Please let me know what you think about this.

Thanking you,

Subarna Sasmal

NYU

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CAGBhpicx%2BUW_-aaqY9nNpuX9B9Xu%3DQ12ESxVeJE%3DwFQpF37pDg%40mail.gmail.com.

[Michele Invernizzi]

Hi,

About LINCS problem: Is your hlda CV based also on the distance between restrained atoms? This could explain the error, because metad would push on this distance, while the LINCS algorithm would try to keep it constant. If this is not the case you could to make the walls more gentle (or the metad, but you have already tried that...)

About the out of grid problem: unfortunately it's enough that one point goes out of range and the simulation crashes. The last printed values in the colvar file can be from a few steps before the crash, at the last flushing (by the way, setting FLUSH STRIDE=1). The only way to fix this is to choose a very broad grid range. There is no need to change the GRID_BIN, 1000 is more than enough even for a grid twice the size (typically a GRID_SPACING of 1/5 the SIGMA is good enough).

Another thing I would suggest is to rescale all the hlda coefficient (as well as the metad SIGMA and GRID) by a common factor of say 1000. This does not change the physics, so it won't alone solve your issues, but might help when setting the metad and wall parameters, as you would be dealing with more natural units.

Hope this helps!

Michele

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CA%2B12EzT-x3pKQu5xcRRRRmWyrPHB8m5_jvxF_o%2BrnudnhYKCBQ%40mail.gmail.com.

[roja rahmani]

Hi,

I have the problem you mentioned (in About LINCS problem). I use time steps 2 fs and the constraints on h-bonds. My peptide has these -h bonds. The distance between the COM of peptide and my nano surface as CV1 and the end to end distance of peptide as CV2. What do you think I should do? Since I use Amber03 FF and dt=2 fs, I have to use the h-bonds for water molecules but how can I differentiate that this constraints should be used only on the h-bonds of water molecules and not my peptides. Do you think I can change the angle lincs-warnangle=30 to more than 30?

Step 42051214, time 84102.4 (ps)  LINCS WARNING in simulation 1

relative constraint deviation after LINCS:

rms 0.000003, max 0.000018 (between atoms 7243 and 7245)

bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length

   7270   7271   32.7    0.1090   0.1090      0.1090

Step 42051222, time 84102.4 (ps)  LINCS WARNING in simulation 1

relative constraint deviation after LINCS:

rms 0.146456, max 1.606307 (between atoms 7270 and 7271)

bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length

   7181   7182   90.0    0.1090   0.1659      0.1090
   7270   7271   90.0    0.1090   0.2841      0.1090
Wrote pdb files with previous and current coordinates

------------------------------------------------------

#RANDOM_EXCHANGES
#RESTART
# slab depth = half slab thickness in nm
depth: CONSTANT VALUES=1.5296195

# Define atom groups
Surface:  GROUP NDX_FILE=index.ndx NDX_GROUP=Surface
Peptide:  GROUP NDX_FILE=index.ndx NDX_GROUP=Peptide
Nterm_CA: GROUP NDX_FILE=index.ndx NDX_GROUP=Nterm-CA
Cterm_CA: GROUP NDX_FILE=index.ndx NDX_GROUP=Cterm-CA

# Make sure that the lipid is not broken by PBC in Gromacs
WHOLEMOLECULES STRIDE=1 ENTITY0=Peptide

# Define virtual atoms and centres-of-mass
Peptide_COM: COM ATOMS=Peptide
Surface_COM: COM ATOMS=Surface
N_CA:        COM ATOMS=Nterm_CA
C_CA:        COM ATOMS=Cterm_CA

# CV1: Distance between Slab COM and Peptide COM (SSD):cv1
# CV2: Peptide end-to end distance (EED):cv2
d: DISTANCE  ATOMS=Surface_COM,Peptide_COM COMPONENTS
d2:  COMBINE ARG=d.z   POWERS=2   COEFFICIENTS=1  PERIODIC=NO
cv1: COMBINE ARG=d2,depth POWERS=0.5,1 COEFFICIENTS=1,-1 PERIODIC=NO
cv2: DISTANCE ATOMS=N_CA,C_CA

METAD ...
   ARG=cv1,cv2 HEIGHT=1.2 SIGMA=0.2,0.2 PACE=500 GRID_MIN=0.0,0.0 GRID_MAX=4.0,3.0 GRID_WSTRIDE=500 GRID_WFILE=GRID LABEL=MW FILE=HILLS
   WALKERS_N=4
   WALKERS_ID=0
   WALKERS_DIR=../HILLS
   WALKERS_RSTRIDE=500
... METAD
wall1: UPPER_WALLS ARG=cv1 AT=3.0 KAPPA=1000000.0 EXP=4 EPS=1.0 OFFSET=0
wall2: UPPER_WALLS ARG=cv2 AT=3.0 KAPPA=1000000.0 EXP=4 EPS=1.0 OFFSET=0

PRINT ARG=cv1,cv2 STRIDE=500 FILE=COLVAR
PRINT ARG=cv1,cv2,MW.bias STRIDE=500 FILE=BIAS
DUMPFORCES STRIDE=500 ARG=cv1,cv2 FILE=FORCES

---------------------------------------------------------------------------------------

Best regards,

Roja