What is required to enabling CMAP scaling for HREX (REST2) with the CHARMM force field

[Chris Neale]

Dear users:

I realize that the charmm CMAP potential is not currently supported for HREX to do REST. However, I wonder if this lack of support means only that the cmap section of the topology files are not processed by "plumed partial_tempering". If this is the case, then I can figure out how to scale the cmap section myself after I have used the plumed script to modify the other parts of the topology. However, if the gromacs mdrun code also needs to be changed to support different cmap potentials among the replicas, then I'll likely pick another force field. I have compelling reasons to use the charmm force field for this particular set of simulations, so I am interested to enable cmap scaling if possible.

Thank you,
Chris.

[Giovanni Bussi]

Hi,

the problem should be only the editing of topologies.

Riccardo Nifosì (in cc) sent me some time ago a script that should correct also cmap. I paste it below, with a small fix provided by Francesco Oteri (also in cc).

awk -v f=$lambda '

BEGIN {cmaptypes=0; pairtypes=0}

{

if (/\[ /) {cmaptypes=0; pairtypes=0};

if (/\[ pairtypes \]/) {pairtypes=1};

if (/\[ cmaptypes \]/) {cmaptypes=1};

if (pairtypes==1) {

sub(/;/, "; ", $0);

print $0;

if ($1!=";" && $1!="[" && NF>0) printf "%s_\t%s_\t%d\t%.12lf  %.12lf    ; scaled \n",$1,$2,$3,$4,$5*f;}

else if (cmaptypes==1) { 

sub(/\\/, "", $0);

if (NF==8) printf "%s %s %s %s %s %s %s %s\\\n",$1,$2,$3,$4,$5,$6,$7,$8;

else if (NF==10) printf "%.9lf %.9lf %.9lf %.9lf %.9lf %.9lf %.9lf %.9lf %.9lf %.9lf\\\n",$1*f,$2*f,$3*f,$4*f,$5*f,$6*f,$7*f,$8*f,$9*f,$10*f; 

else if (NF==6) printf "%.9lf %.9lf %.9lf %.9lf %.9lf %.9lf\n",$1*f,$2*f,$3*f,$4*f,$5*f,$6*f;

else print $0;} 

else {print $0;}

}

'

Remember to check that energy terms are properly scaled (e.g. with a rerun).

Concerning the modification to mdrun, they should be compatible with cmaps. Anyway, be always careful with this patched gromacs, it is not as robust as normal gromacs+plumed.

Regards,

Giovanni

--
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 http://groups.google.com/group/plumed-users.
To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/154ee874-d76f-4453-a957-acea59d225b2%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Chris Neale]

Awesome! Thank you very much. I'll be sure to test things thoroughly.

[Shyno Mathew]

Hi all,

I am trying to fix the cmap scaling

this fix goes to the partial_tempering.sh ?

thanks,

Shyno

[Riccardo NIFOSì]

Hi,

the awk script to perturb the cmap terms goes as a post-processing after partial_tempering.sh.

so, for example 

plumed partial_tempering $lambda < topol.top | 

awk -v f=$lambda ' BEGIN {cmaptypes=0; pairtypes=0}

{

if (/\[ /) {cmaptypes=0; pairtypes=0};

if (/\[ pairtypes \]/) {pairtypes=1};

if (/\[ cmaptypes \]/) {cmaptypes=1};

if (pairtypes==1) {

print $0;

if (substr($1,0,1)!=";" && substr($1,0,1)!="[" && NF>0) printf "%s_\t%s_\t%d\t%.12lf  %.12lf    ; scaled \n",$1,$2,$3,$4,$5*f;}

else if (cmaptypes==1 && substr($1,0,1)!=";") { 

sub(/\\/, "", $0);

if (NF==8) printf "%s %s %s %s %s %s %s %s\\\n",$1,$2,$3,$4,$5,$6,$7,$8;

else if (NF==10) printf "%.9lf %.9lf %.9lf %.9lf %.9lf %.9lf %.9lf %.9lf %.9lf %.9lf\\\n",$1*f,$2*f,$3*f,$4*f,$5*f,$6*f,$7*f,$8*f,$9*f,$10*f; 

else if (NF==6) printf "%.9lf %.9lf %.9lf %.9lf %.9lf %.9lf\n",$1*f,$2*f,$3*f,$4*f,$5*f,$6*f;

else print $0;} 

else {print $0;}

}

'  > topol$i.top  

###################

(it is slightly modified with respect to the previous version, to fix a bug) 

This perturbs cmap terms for the WHOLE protein part (this is because the cmap types are linked to the bond types of the involved atoms, which are not perturbed in the original plumed hrex scheme).

So in case you need to perturb only a subregion of your protein, this is not going to work. 

Regards,

Riccardo

[Shyno Mathew]

Thanks Riccardo.

I am perturbing the whole protein, so this works. Now for lambda=0.5, I am getting CMAP dih. values as half compared to those for lambda=1

thanks,

Shyno

--

Sincerely
Shyno Mathew

[Chris Neale]

Thank you Riccardo.

I'll note that if you have a system that includes [ pairtypes ] and you don't use the script above to modify the output from the plumed partial_tempering run, then you'll wind up in a situation where if you script everything and your lambda=1 replica was actually obtained by running "plumed partial_tempering 1 < topol.top > topol0.top" then you don't get the right force field. I presume that the same goes for anything including cmap. The issue is that the hot part gets "_" added to the type and so it doesn't get recognized in places where the [ pairtypes ] and [ cmap ] have not been updated -- and the result is that these interactions get completely turned off. So it's not at all just a question of "do I want to scale the 1-4 LJ interactions or the cmap?" This is serious enough that I'd almost call it a bug in the plumed/gromacs/hrex concept and even though it might technically not be a bug in the hrex script I'd suggest that the hrex script incorporate this.

In any event, I have just now verified that for a complex system the hrex 2.2.2 -hrex version gives identical energies at lambda=1 and with an unperturbed topology that never touched plumed. However, since my system had [ pairtypes ], things were not correct until I used the above script from Riccardo. 

Thank you,

Chris.

[Giovanni Bussi]

Hi all,

If someone sends me an example (topol.to, conf.gro and mdrun.mdp) I can incorporate the changes that Riccardo made in the distribute script to make life easier.

Giovanni

Inviato da iPhone

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/2532a3a0-3baf-4c6b-9bff-0c3b62fec16c%40googlegroups.com.

[Chris Neale]

Giocanni: I'll send you a set of files now

Riccardo: can you please explain why you don't need to add underscores in the cmap part? 

Thank you,

Chris.

[Riccardo NIFOSì]

Dear Chris,

my choice was one of convenience. In GROMACS the non-bonded and bonded atom types are treated differently. In the example I give below,  the first two columns the [atomtypes] section are "name" and "bond_type". The  first is the non-bonded type, and the second the bonded one. These are used for default assignment of bonds and angle interactions. In Giovanni's implementation bond and angle interactions are not perturbed, and indeed the bond type column is unchanged (so that for example perturbed N3_ has the same bond type, N3). The CMAP interactions are described by bonded types, so it was easier for the kind of application I had in mind to just scale all the CMAP terms instead of defining new bondtypes (and consequently re-listing all angle and bond terms for perturbed types and also all possible cross terms, which should anyway be easy to do). 

Hope this helps!

Riccardo

[ atomtypes ] 

;name   bond_type     mass     charge   ptype   sigma         epsilon       Amb

N3 N3 0.00000 0.00000 A 3.25000e-01 7.11280e-01 ; 1.82  0.1700

N3_ N3 0.00000 0.00000 A 3.25000e-01 0.699079  ; scaled

H H 0.00000 0.00000 A 1.06908e-01 6.56888e-02 ; 0.60  0.0157

H_ H 0.00000 0.00000 A 1.06908e-01 0.064562  ; scaled

[Chris Neale]

Ah, great. I didn't realize this. Thank you Riccardo.

[Chris Neale]

Dear all:

I think this treatment of pairtypes is still going to miss cases in which one member of the pair is in the hot region.

To fix that, I suggest to replace the line:

      if (substr($1,0,1)!=";" && substr($1,0,1)!="[" && NF>0) printf "%s_\t%s_\t%d\t%.12lf  %.12lf    ; scaled \n",$1,$2,$3,$4,$5*f;}

with:

      

      if (substr($1,0,1)!=";" && substr($1,0,1)!="[" && NF>0){

         printf "%s_\t%s_\t%d\t%.12lf  %.12lf    ; scaled \n",$1,$2,$3,$4,$5*f;

         printf "%s_\t%s\t%d\t%.12lf  %.12lf    ; scaled \n",$1,$2,$3,$4,$5*sqrt(f);

         printf "%s\t%s_\t%d\t%.12lf  %.12lf    ; scaled \n",$1,$2,$3,$4,$5*sqrt(f);

      }

I realize that you designed this script to decouple an entire protein (since otherwise there are issues with cmap), but I think that the above change still generalizes the script in a useful way.

Chris.

[Giovanni Bussi]

Hi all,

Chris: you are right. I am integrating the changes in the partial_tempering script and did it in the way you suggest.

I am trying to make the cmap stuff but I have some doubt. The script by Riccardo makes some assumption about how many fields per line are present. I am not sure this is a rule in the topol.top syntax. I will check the way topol.top is read in gromacs.

Giovanni

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/3f1d03bf-4ad4-41e1-b6ff-678daf9f3c7e%40googlegroups.com.

[rajeswar...@gmail.com]

Hello all,
I am intending to do REST simulation in gromacs with charmm36m forcefield. I am using plumed 2.3.2 version. While generating the scaled topologies using partial tempering, do I still need to post process the output using the awk script provided here? I guess the [pairtypes] have already been scaled by partial tempering. Hence my question is if I only correct the cmap potential on the output of partial tempering, is that sufficient? Please suggest me.

Thank you,
Rajeswari A

[Chris Neale]

It's been a while, but I think the answer is yes. The script is required to do the pairtypes and dihedrals as well. You can test by running the script and then doing a diff on the input and output to see what changed. (Or meld might be easier to use ... it's a visual version of diff)

--
You received this message because you are subscribed to a topic in the Google Groups "PLUMED users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/plumed-users/Txlf8Icwr2k/unsubscribe.
To unsubscribe from this group and all its topics, send an email to plumed-users+unsubscribe@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/7dcb61f9-c199-420a-9534-973466de96df%40googlegroups.com.

[Rajeswari A.]

Thank you Chris for your prompt reply. The output from partial tempering I get the pairtypes section already scaled. For your reference, here is my output from partial tempering. Only the pairtypes section is given.
[ pairtypes ]
CP1 C 1 0.347450500075 0.138767581228
CP1_ C 1 0.347451 0.148702  ; scaled
CP1 C_ 1 0.347451 0.148702  ; scaled
CP1_ C_ 1 0.347451 0.159348  ; scaled
CP1 CA 1 0.346773417049 0.110698234855
CP1_ CA 1 0.346773 0.118623  ; scaled
CP1 CA_ 1 0.346773 0.118623  ; scaled
CP1_ CA_ 1 0.346773 0.127116  ; scaled

The corresponding commands in the partial tempering script can also be find here,

# ATOMTYPES (PAIRS)
  } else if((rec=="pairtypes" && NF>=5) || rec=="nonbond_params" && NF>=5){
    scale2=1.0; # scaling for second to last column
    if(combrule==1) scale2=scale;
    print $1,$2,$3,$4,$5,comments
    print $1""suffix,$2,$3,sqrt(scale2)*$4,sqrt(scale)*$5," ; scaled";
    print $1,$2""suffix,$3,sqrt(scale2)*$4,sqrt(scale)*$5," ; scaled";
    print $1""suffix,$2""suffix,$3,scale2*$4,scale*$5," ; scaled";

On top of this, if I use the awk command provided by Riccardo, it is scaling also the already scaled pairtypes. Example output is here.
CP1 C 1 0.347450500075 0.138767581228
CP1_    C_      1       0.347450500075  0.071696629890    ; scaled
CP1_    C       1       0.347450500075  0.099745515749    ; scaled
CP1     C_      1       0.347450500075  0.099745515749    ; scaled
CP1_ C 1 0.347451 0.0997455  ; scaled
CP1__   C_      1       0.347451000000  0.051535208248    ; scaled
CP1__   C       1       0.347451000000  0.071696618570    ; scaled
CP1_    C_      1       0.347451000000  0.071696618570    ; scaled
CP1 C_ 1 0.347451 0.0997455  ; scaled
CP1_    C__     1       0.347451000000  0.051535208248    ; scaled
CP1_    C_      1       0.347451000000  0.071696618570    ; scaled
CP1     C__     1       0.347451000000  0.071696618570    ; scaled


Does it needed to do that way only? Or can i comment out the pairtype scaling from the awk script provided by Riccardo?

Thanks once again for your help!,

Rajeswari A,

SERB National Postdoctoral Fellow

Molecular Biophysics,

Indian Institute of Science,
Bangalore, India - 560 012

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CAOge3ziCnt%2BZ1%2BStjFChxKug8CJArZqTgUkEk5fEqC6YndOeGA%40mail.gmail.com.

[davtya...@gmail.com]

I am trying to run REST2 simulations using this thread. So far, based on energy differences between unscaled simulation and lambda=0.5 case, it seems to me that the pairtypes modification is no longer necessary and only the CMAP part of the script needs to be applied. Can anyone confirm this? 

I do run into a problem when trying to run the REST2. The simulations seem to proceed normally till the first exchange (or if no attempt for exchange is made). However, I get the following error as soon the first attempt is made:

step 500: One or more water molecules can not be settled.

Check for bad contacts and/or reduce the timestep if appropriate.

At the same time, the simulations seems to run normally for any of the modified systems when -plumed flag is not used (i.e. when I run a regular MD simulation with gromacs using the scaled topology files).

I use Gromacs/2018.4 + Plumed2.5.0. Have anyone encountered similar problems? 

To unsubscribe from this group and all its topics, 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/7dcb61f9-c199-420a-9534-973466de96df%40googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to a topic in the Google Groups "PLUMED users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/plumed-users/Txlf8Icwr2k/unsubscribe.

To unsubscribe from this group and all its topics, 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.