Restraint Bias in GROMACS Steepest Descent

[Kelvin Azevedo Santos]

Hi all,

I would like to run mdrun with steep integrator using plumed bias. I have found here (https://groups.google.com/forum/#!msg/plumed-users/0vu_lmUdiq0/TyCHVnV1MAAJ) that it used not to be supported, but that after version 2.1 it was implemented. However, I could not find out how to make it work.

I run "gmx mdrun -s minim.tpr -plumed plumed.dat" without any warnings or errors, but the result is the same as it would have without adding the plumed file, just as it was reported in the other thread to be before version 2.1 implemented this functionality. What is it that could be wrong?

My gromacs is version 2016.5

My plumed is version 2.4

Thank you for your help,

Kelvin A. Santos

[Carlo Camilloni]

More details are needed, as said it is implemented since version 2.1

Sent from my iPhone

--
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 post to this group, send email to plumed...@googlegroups.com.
Visit this group at https://groups.google.com/group/plumed-users.
To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/2b9bb118-1b59-4f74-8a0f-2f19bc0d6b10%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Kelvin Azevedo Santos]

Thank you for your quick feedback. Let me give you as much information as I can. First I have a reference pdb, called "modes.pdb", containing the eigenvalues for a PCA decomposition in the first 8 modes (therefore 9 frames). This is the pdb that will be used as reference for the PCAVARS command in our .dat. Out "modes.pdb" has about 10000 atoms, then "END", "REMARK TYPE=OPTIMAL", then another round of 10000 atoms, and so on, until the last "END".

Out plumed.dat file is this one:

PCAVARS REFERENCE=modes.pdb LABEL=pca

RESTRAINT ARG=pca.eig-1 KAPPA=100.0 AT=+1.0 LABEL=lr

PRINT ARG=lr.bias FILE=result.data STRIDE=1

We then have this minim.mdp:

integrator = steep

emtol = 100.0

emstep = 0.001

nsteps = 50000

nstlist = 10

cutoff-scheme = Verlet

ns_type = grid

coulombtype = PME

rcoulomb = 1.4

rvdw = 1.4

pbc = xyz

We then have the previously minimized "protein.gro" file, and the topology (topol.top). I run the following commands:

gmx grompp -f minim.mdp -c protein.gro -p topol.top -o minim.tpr

gmx mdrun -s minim.tpr -o traj.trr -c mode1.gro -plumed plumed.dat

It then runs for 12 steps and stops, because the protein.gro file is already minimized (that is, already at=0 for every eigenvector). Of course, it would have run for many more steps if it really had the bias working correctly. The output is exactly the same as the one generated by not adding the "-plumed plumed.dat".

In the md.log file, I get these messages from PLUMED:

PLUMED: PLUMED is starting

PLUMED: Version: 2.4.2 (git: Unknown) compiled on Aug 21 2018 at 21:13:27

(...)

PLUMED: FILE: plumed.dat
PLUMED: Action PCAVARS
PLUMED:   with label pca
PLUMED:   found 8 eigenvectors in file modes.pdb
PLUMED:   added component to this action:  pca.eig-1
PLUMED:   added component to this action:  pca.eig-2
PLUMED:   added component to this action:  pca.eig-3
PLUMED:   added component to this action:  pca.eig-4
PLUMED:   added component to this action:  pca.eig-5
PLUMED:   added component to this action:  pca.eig-6
PLUMED:   added component to this action:  pca.eig-7
PLUMED:   added component to this action:  pca.eig-8
PLUMED:   added component to this action:  pca.residual
PLUMED: Action RESTRAINT
PLUMED:   with label lr
PLUMED:   with arguments pca.eig-1
PLUMED:   added component to this action:  lr.bias
PLUMED:   at 1.000000
PLUMED:   with harmonic force constant 100.000000
PLUMED:   and linear force constant 0.000000
PLUMED:   added component to this action:  lr.force2
PLUMED: Action PRINT
PLUMED:   with label @2
PLUMED:   with stride 1
PLUMED:   with arguments lr.bias
PLUMED:   on file result.data
PLUMED:   with format  %f
PLUMED: END FILE: plumed.dat
PLUMED: Timestep: 0.001000
PLUMED: KbT has not been set by the MD engine
PLUMED: It should be set by hand where needed
PLUMED: Relevant bibliography:
PLUMED:   [1] Tribello, Bonomi, Branduardi, Camilloni, and Bussi, Comput. Phys. Commun. 185, 604 (2014)
PLUMED: Please read and cite where appropriate!
PLUMED: Finished setup
Started Steepest Descents on rank 0 Fri Oct 26 10:49:44 2018

(...)

In the end I have more lines from PLUMED:

Finished mdrun on rank 0 Fri Oct 26 10:51:41 2018
PLUMED:                                               Cycles        Total      Average      Minumum      Maximum
PLUMED:                                                    1   117.475397   117.475397   117.475397   117.475397
PLUMED: 1 Prepare dependencies                            13     0.000039     0.000003     0.000001     0.000022
PLUMED: 2 Sharing data                                    13     0.001295     0.000100     0.000040     0.000805
PLUMED: 3 Waiting for data                                13     0.000021     0.000002     0.000001     0.000003
PLUMED: 4 Calculating (forward loop)                      13   117.026780     9.002060     8.871689     9.228949
PLUMED: 5 Applying (backward loop)                        13     0.001582     0.000122     0.000116     0.000140
PLUMED: 6 Update                                          13     0.000163     0.000013     0.000009     0.000048

Here is the result.data:

#! FIELDS time lr.bias
 0.000000 50.000000
 0.001000 50.000000
 0.002000 50.000000
 0.003000 50.000000
 0.004000 50.000000
 0.005000 50.000000
 0.006000 50.000000
 0.007000 50.000000
 0.008000 50.000000
 0.009000 50.000000
 0.010000 50.000000
 0.011000 50.000000
 0.012000 50.000000

I am sorry if this information is not enough, or is irrelevant, or both. I am quite a begginer with Plumed yet, but I am learning more every day. Again thank you for your feedback and please tell me, if the above information is not enough for starting a diagnosis, how can I generate the required information so that you can help me solve this problem. Thank you,

Kelvin Santos

[Carlo Camilloni]

Or 12 steps are enough to reach the minimization target of a maximum force of less than 100. I am not sure about the pca because Is something I don’t use, but as a Simple test increase the kappa to 1000000 or so and repeat the minimization 

C

Sent from my iPhone

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/11568ac0-7a6e-4991-a2c3-4ae8988fecba%40googlegroups.com.

[Kelvin Azevedo Santos]

I have repeated the tests with KAPPA=1000000 and nothing changed. The result.data log now shows that "lr.bias" is 500000.0, and the "md.log" file shows the new kappa, but otherwise we get no results. If I change the "AT=" parameter, or try a linear restraint, that changes nothing. If instead of using PCAVARS I use a normal, distance between atoms restraint, again I get no result, so the problem is nothing related to PCAVARS.

I have tried as well to decrease the maximum force to 1, and in this case Steepest Descent converges to machine precision within 39 steps, with no other difference. The same happens when I change other parameters of the mdp file.

I would be very much grateful if you could tell me of a test which would allow me to progress past this point... should I check the "minimize.c" which was patched? Is it so much difficult for a plumed begginer (with reasonable C-language skills) to understand where the problem is? I wish I could do some debugging.

[Carlo Camilloni]

What I find weird is that the value of the CV is not changing. Is the same happening also in the case of a distance?

You can debug the code for sure. 

The patch for the minimisation is the minimise.cpp file in mdlib folder and each part of the patch is between a 

/* PLUMED */

...

/* END PLUMED */ 

comment.

C

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/32334bf0-1447-4b21-a1f7-dee61521aa3a%40googlegroups.com.