hbonds as colvar

[Erica Prates]

Dear plumed users,

I want to compute the free energy surface across the radius of gyration and intramolecular hydrogen bonds of a polypeptide chain. However, I am not sure about the best way to count the hydrogen bonds in plumed 2.0. I am trying to use COORDINATION, with a gaussian switching function, considering its peak at 0.2 nm and a small width of 0.05 nm. But in COLVAR file, the resulting number of hbonds is too high, suggesting that even the covalent bonds involving the hydrogens are being counted. So, is there another better way I could do this hbond calculation? My input looks like this now (groupA refers to the hydrogens and groupB, to the hbond acceptors):

WHOLEMOLECULES ENTITY0=9,23,37,61,67,81,103,109,123,145,151,167,181,203,209,223,235,247,258
rg: GYRATION TYPE=RADIUS ATOMS=9,23,37,61,67,81,103,109,123,145,151,167,181,203,209,223,235,247,258
c: COORDINATION GROUPA=8,18,22,32,36,66,76,80,90,108,118,122,132,150,166,176,180,190,208,218,222,234,246,253,257,266,267,271  GROUPB=6,7,17,20,21,31,34,35,50,51,64,65,75,78,79,89,92,93,106,107,117,120,121,131,134,135,148,149,164,165,175,178,179,189,192,193,206,207,217,220,221,229,230,232,233,241,242,244,245,252,255,256,264,265,269,270 SWITCH={GAUSSIAN D_0=0.2 R_0=0.05}

metad: METAD ARG=* PACE=500 HEIGHT=1.2 SIGMA=0.05,0.05 FILE=HILLS-RG-HB BIASFACTOR=10.0 TEMP=298

PRINT STRIDE=100 ARG=rg,c,metad.bias FILE=COLVAR-RG-HB

[Gareth Tribello]

Hello

If you are willing to change to the master branch of the code you could try the command HBOND_MATRIX:

http://plumed.github.io/doc-master/user-doc/html/HBOND_MATRIX.html

For this you specify switching functions on the OO distance the OH distance and the OHO angle.  To get something equivalent to the COORDINATION (but with hydrogen bonds) you will need a set of commands something like this:

mat: HBOND_MATRIX DONORS ACCEPTORS HYDROGENS

sums: ROWSUMS MATRIX=max SUM

PRINT ARG=sums.sum

I hope this helps

Gareth 

--
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/966a6080-0c75-4e1e-9361-bf5abe262d3c%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Erica Prates]

Hello Gareth

Thank you very much for your time and help.

So, I uploaded plumed at https://github.com/plumed/plumed2/tree/40327addd5220fe1b6d6e9e0f16521fe93d745b7  . Is this what you meant with "change to the master branch"?
Then, I configured that with -enable-modules=adjmat . That worked perfectly.
However, when I tried to  patch this version to gromacs, it claimed the plumed wasn't configured with MPI. That is quite strange, since, yes, it was configured with MPI... I configured plumed exactly the same way as previously.  Have you faced this same issue?


Best
Erica

[Gareth Tribello]

Hello

I am not sure if what you have downloaded is what I meant.  Can you try the following command:

git clone https://github.com/plumed/plumed2.git

Then do configure in the way you did before and make plumed.  

I am not sure why you get MPI issues at the moment.  I guess that this is during runtime?  It is indeed strange and I haven’t seen it before.  Try the above though and see what happens.  

Gareth

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/910d5e42-eec2-4d4e-941c-e8ea797d141a%40googlegroups.com.

[Erica Prates]

Hello Gareth,


Thank you very much for the help! The installation worked properly after using that command you suggested.
However, how to write correctly an input for HBOND_MATRIX is not clear at all to me yet. When I tried following the example you gave me, I get the error:


message: ERROR in input to action ROWSUMS with label sums : max does not calculate an adjacency matrix


I think I can set the action SUM within the mat descriptors. So, I wrote the input as:


WHOLEMOLECULES ENTITY0=1-275
rg: GYRATION TYPE=RADIUS  ATOMS=... 
mat: HBOND_MATRIX SUM HYDROGENS=... ACCEPTORS=... DONORS=... SWITCH={RATIONAL R_0=0.3}  ASWITCH={RATIONAL R_0=30} HSWITCH={RATIONAL R_0=0.105} 

metad: METAD ARG=* PACE=500 HEIGHT=1.2 SIGMA=0.05,0.05 FILE=HILLS-RG-HB BIASFACTOR=10.0 TEMP=298

PRINT STRIDE=100 ARG=rg,mat.sum,metad.bias FILE=COLVAR-RG-HB


It started running, and the description in the log file seems to be right. However, in the first step of the simulation, the value computed for mat.sum is not recognized (in COLVAR file it gives -nan), and the simulation stops. Do you know what can be wrong in my input file?
In the user guide of HBOND_MATRIX, it is said we should specify ATOMS OR ACCEPTORS DONORS OR HYDROGENS. It is confused.. Also, do you know if it is right the angle in ASWITCH be specified in degrees?
I pasted bellow the log lines.


Best,
Erica


Log:

PLUMED: Action HBOND_MATRIX
PLUMED:   with label mat
PLUMED:   involving atoms 7 17 21 31 35 65 75 79 89 107 117 121 131 149 165 175 179 189 207 217 221 233 245 252 256 265 270
PLUMED:   involving atoms 6 17 20 31 34 50 64 75 78 89 92 106 117 120 131 134 148 164 175 178 189 192 206 217 220 229 230 232 241 242 244 252 255 264 269
PLUMED:   involving atoms 8 18 22 32 36 66 76 80 90 108 118 122 132 150 166 176 180 190 208 218 222 234 246 253 257 266 267 271
PLUMED:   involving hydrogen atoms : 8 18 22 32 36 66 76 80 90 108 118 122 132 150 166 176 180 190 208 218 222 234 246 253 257 266 267 271
PLUMED:   atoms of type 1 and 1 must be within 0.3.  Using rational swiching function with parameters d0=0 nn=6 mm=12
PLUMED:   for atoms of type 1 and 1 the OH distance must be less than 0.105.  Using rational swiching function with parameters d0=0 nn=6 mm=12
PLUMED:   for atoms of type 1 and 1 the OOH angle must be less than 30.  Using rational swiching function with parameters d0=0 nn=6 mm=12
PLUMED:   added component to this action:  mat.sum
PLUMED:   value mat.sum contains the sum of all the values
PLUMED: Action METAD
PLUMED:   with label metad
PLUMED:   with stride 1
PLUMED:   with arguments rg mat.sum
PLUMED:   Gaussian width  0.050000 0.050000  Gaussian height 1.200000
PLUMED:   Gaussian deposition pace 500
PLUMED:   Gaussian file HILLS-RG-HB
PLUMED:   Well-Tempered Bias Factor 10.000000
PLUMED:   Hills relaxation time (tau) 18.582824
PLUMED:   KbT 2.477710
PLUMED:   added component to this action:  metad.bias
PLUMED:   added component to this action:  metad.work
PLUMED:   Bibliography [2][3]
PLUMED: Action PRINT
PLUMED:   with label @4
PLUMED:   with stride 100
PLUMED:   with arguments rg mat.sum metad.bias
PLUMED:   on file COLVAR-RG-HB
PLUMED:   with format  %f
PLUMED: END FILE: plumed.dat
PLUMED: Timestep: 0.002000
PLUMED: KbT: 2.478971
PLUMED: Relevant bibliography:
PLUMED:   [1] Tribello, Bonomi, Branduardi, Camilloni, and Bussi, Comput. Phys. Commun. 185, 604 (2014)
PLUMED:   [2] Laio and Parrinello, PNAS 99, 12562 (2002)
PLUMED:   [3] Barducci, Bussi, and Parrinello, Phys. Rev. Lett. 100, 020603 (2008)
PLUMED: Please read and cite where appropriate!
PLUMED: Finished setup
Started mdrun on rank 0 Tue Mar  8 08:11:17 2016
           Step           Time         Lambda
              0        0.00000        0.00000

   Energies (kJ/mol)
           Bond          Angle    Proper Dih.  Improper Dih.          LJ-14
    2.09437e+02    6.39163e+02    8.56944e+02    3.31116e+01    2.75620e+02
     Coulomb-14        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential
    2.76741e+03    8.21471e+04   -5.95924e+05    3.42451e+03   -5.05571e+05
    Kinetic En.   Total Energy    Temperature Pressure (bar)   Constr. rmsd
           -nan           -nan           -nan           -nan            nan

[Gareth Tribello]

Hello Erica

In that input you sent if you did the command

mat: HBOND_MATRIX SUM HYDROGENS=... ACCEPTORS=... DONORS=... SWITCH={RATIONAL R_0=0.3}  ASWITCH={RATIONAL R_0=30} HSWITCH={RATIONAL R_0=0.105}  

ROWSUMS MATRIX=max

Then that won’t work because the label of HBOND_MATRIX is mat.  

Other things is that for ASWITCH, which is the switching function on the angle, the units are radians and not degrees - remember though that you can put numbers in terms of pi in your plumed input.  So 0.1667pi is OK for 30 degrees.  

I think you are right though that the SUM should also work in this case.  

If it doesn’t work with the change of units from degrees to radians and you want me to take a look send me the input you have and some frames and I will see what I can do.

Thanks

Gareth

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/7e27c56b-e353-48e4-a7ab-34239339fd92%40googlegroups.com.

[Erica Prates]

I sent my files to Gareth in private, so I paste here Gareth's e-mail that solved my end question. This might be helpful to someone else.

I had a look at this and there is a problem in the input file.  First of all the idea of the CV is as follows.  A hydrogen bond looks something like this:

D-H——A

You thus need three coordinates to describe the geometry of the hydrogen bond.  For the coordinate you state that to have a hydrogen bond the distance between the atom D and the atom A must be in a certain rate.  The distance between the atom D and the atom H must be in a certain range and the DHA angle must be in a certain range.  This is what SWITCH, HSWITCH and the ASWITCH commands in your input are doing.  

The atoms like the D in the above are the set of DONOR atoms in your molecule.  In other words, the donor atoms have a hydrogen covalently bonded to them.  The As meanwhile are the atoms in the ACCEPTOR group in the plumed input file.  To operate in the way you are operating there should be no atoms that are simultaneously in the DONOR and ACCEPTOR groups.  These two groups of atoms should be mutually exclusive.  Either an atom has a hydrogen covalently bound to it or it does not.

The issue in your case is that the two groups are not mutually exclusive.  This is causing all kinds of problems because there are places in the code where you are dividing by zero, which is the reason you are getting nan.

You can get round this problem by operating using the ATOMS keyword instead of DONOR and ACCEPTOR.  In this case both the acceptor and DONOR groups are assumed to contain the exact same set of atoms. 

[thidar...@gmail.com]

Hi Erica and Gareth

I am going to use Hbond_matrix as one of CVs for my protein like you did.

Could you please share your plumed.dat which include the command for HBOND_MATRIX that worked for your system?

I am really struggling..

Regards,

Tar

[Gareth Tribello]

Hello Tar

Here is an example input for HBOND_MATRIX:

mat: HBOND_MATRIX ATOMS=1-192:3 HYDROGENS=2-192:3,3-192:3 SWITCH={RATIONAL R_0=3.20} HSWITCH={RATIONAL R_0=2.30} ASWITCH={RATIONAL R_0=0.167pi} SUM

ROWSUMS MATRIX=mat MEAN LABEL=rsums

COLUMNSUMS MATRIX=mat MEAN LABEL=csums

PRINT ARG=rsums.mean,csums.mean FILE=colvar

So this is calculating the number of hydrogen bonds formed by each of the atoms from 1-192:3 and the hydrogens in those lists.  It is a system of water molecules.  If you want to look more closely all the files for this example are in the plumed directory you downloaded in 

plumed2/regtest/adjmat/rt-hbond

If you go that directory and type make the calculation will run automatically and you can look at the input and output.

I hope this helps

Gareth

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/c8268003-343b-4784-91a0-e20a74b449ba%40googlegroups.com.

[Thidarat Wongpinyochit]

Dear Gareth 

Thank you so much for sharing! That's very helpful!

However, I have question about the ATOMS and HYDROGENS. They includes all atom in the system? If yes, will it make a lot of hydrogen bond like including bonds in the same residue?

In my cases, I have 304 residues (300 amino acids + 4 terminal cuz it contains 2 chains) which contain 3194 atoms.I would like to determine the number of h-bonding between O of residue i and N of residues i+4 for example..

I will use radius of gyration and the number of h-bond as CVs in metadynamics.

Could you please clarify this..

Thank you in advance!!!

Regards,

tar

To unsubscribe from this group and stop receiving emails from it, send an email to plumed-users+unsubscribe@googlegroups.com.
To post to this group, send email to plumed-users@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/c8268003-343b-4784-91a0-e20a74b449ba%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/gZRPu_8HX7s/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/394491F7-172F-4CC5-87D4-5E10DA03369F%40gmail.com.

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

--

Thidarat Wongpinyochit

PhD student 

Strathclyde Institute of Pharmacy and Biomedical Sciences

University of Strathclyde

Tel: +447928348679

[Gareth Tribello]

Hello

The example I sent you uses all the oxygen atoms in the water molecules in the system yes.  It doesn’t have to be so expensive though as plumed will use various tricks to calculate things.  If you know what hydrogen bonds are going to form you can use DONORS and ACCEPTORS.  This will tell you if any of the DONOR atoms donates a hydrogen bond to the ACCEPTOR atoms.  You need to work out what atoms are going to potentially form hydrogen bonds though as you know your system and what it is you want to do.

Gareth

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/CABFDis%2BuZ2SDWqQHoD6X%3DzNaG4vU03%3DUQ1H1LgLYFJWjmCnehw%40mail.gmail.com.

[Thidarat Wongpinyochit]

Dear Gareth

Thank you so much!

To unsubscribe from this group and stop receiving emails from it, 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/CABFDis%2BuZ2SDWqQHoD6X%3DzNaG4vU03%3DUQ1H1LgLYFJWjmCnehw%40mail.gmail.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/gZRPu_8HX7s/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/BA851FB0-D5B9-419B-947A-4CD42D19C7CB%40gmail.com.

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