about using ENERGY as CV in LAMMPS

[James Yuan]

Hi All:

I'm new to the PLUMED and I'm trying to get ENERGY output as CV on the fly in LAMMPS (without any bias). It gives me a list of errors. However, when I switch to DISTANCE as CV, I'm able to output it smoothly. I wonder if PLUMED in LAMMPS supports ENERGY as CV for this case? If not, is there any way to code it? 

BTW: my lammps version is 1Feb14 and plumed version is 2.1.0. 

Thanks! 

James

[Giovanni Bussi]

Hi,

ENERGY is currently not supported in LAMMPS. I am not a LAMMPS expert, and I am not sure if it is possible to implement it. If you know how to access to the total energy in a LAMMPS fix, we can tell you how to pass it to 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/dce7214f-71b2-4cb6-9bb9-619c04623967%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[James Yuan]

Hi Giovanni:

Thanks for you quick reply. In lammps, we can quickly invoke a compute to compute the total energy by referring to pair style. The beauty lies in its generality, i.e. a simple code can deal with all compatible potentials in lammps. I can easily incorporate this into "fix plumed" if you can show me how CV variables are defined and used by Plumed.

Best,
James

[Giovanni Bussi]

Hi,

to pass the total energy to plumed you should just add am instruction in the fix_plumed.cpp file.

After line:

p->cmd("setVirial",&virial[0][0]);

add:

p->cmd("setEnergy",&energy);

where energy is a double precision variable storing the energy.

In case to fill the energy variable you have to do something expensive, you can ask PLUMED is energy is actually required:

int needsEnergy;

p->cmd("isEnergyNeeded",&needsEnergy);

if(needsEnergy){

// do expensive call

  p->cmd("setEnergy",&energy);

}

To validate the implementation you can try to add a bias on the ENERGY and see if its distribution is properly affected.

E.g., compute F(E)=-kT*log P(E), where P(E) is the histogram of the energy.

Then add a bias:

e: ENERGY

RESTRAINT AT=0 ARG=e KAPPA=0 SLOPE=1

and recompute the new F'(E). Difference between F' and F should be a straight line with slope 1.

Let us know if this change works. If so we can add it to the official version.

Thanks!

Giovanni

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/09265ee6-6726-4b00-8757-078dda84c707%40googlegroups.com.

[James Yuan]

Hi Giovanni:

I've successfully added in the code as you suggest and my current implementation can set energy value to that from LAMMPS' output. However, I found that when I did a simple on-the-fly calculation by a "fix plumed" command, the energy actually is two times the correct value after I've taking into consideration of unit conversion. So I'm wondering what kind of energy we need to set by "p->cmd("setEnergy",&energy);" or alternatively what kind of "ENERGY" we output by the following plumed.dat?  BTW, I currently include only potential energy but kinetic energy can be easily added in. 

My plumed.dat for "fix plumed" is attached below:

----

pe: ENERGY 

PRINT ARG=pe STRIDE=1 FILE=PE.COLVAR

-----

Best,

James

[Giovanni Bussi]

Hi,

which units are you using in LAMMPS? Conversion factors should be set also in fix_plumed.cpp. Look for "LAMMPS units" and below.

Giovanni

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/e2f20e4e-6194-4228-99d6-152ed3248120%40googlegroups.com.

[James Yuan]

Hi Giovanni:

I've used metal units and I've successfully got the correct output for potential energy. However, I found that I could get the correct output when I used one core to the lammps (with plumed patched). Whenever, I use N core to run the same input script, the plumed gives me N times the correct value. Is this a potential bug of plumed? Please take your time to check it out. I've attached my input for your convenience (a FCC Cu metal with EAM potential). 

Best,

James

[Giovanni Bussi]

Hi,

PLUMED assumes that every processor stores part of the energy and the sum is the total energy. This is something that originates from the implementation in GROMACS.

To fix this, just divide the energy times by number of processors (or set it to zero on non-root processors) before passing it to PLUMED.

Giovanni

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/8bc090f2-00f5-4950-b6a3-54c78333056f%40googlegroups.com.

[James Yuan]

Hi Giovanni:

Thanks for pointing this out. Now I'm confident that energy as CV is working in LAMMPS. I'll contribute this part to future update of Plumed. 

Best,

James

[davide.mi...@gmail.com]

Dear Giovanni and James,

I know this thread has been untouched for a while but I was wondering if you could point me on the right direction.
I'd like to implement the potential energy as CV myself.
I have added the piece of code suggested by Giovanni but I can't figure out a simple way to access the potential energy of the system from fix_plumed.cpp
so that it can be assigned to a variable energy.

I guess this is more a LAMMPS question rather than a PLUMED question.

Any suggestion would be very useful.

Thanks in advance.

Davide

[João Fernandes]

Hello guys.

I made as suggested changes, a newer error.

terminate called after throwing an instance of 'PLMD::Exception'

  what():  

+++ Internal PLUMED error

+++ file Atoms.h, line 197

+++ message: assertion failed collectEnergy && energyHasBeenSet

Any suggestion?

I am using plumed2.2 in VMD

DISTANCE ATOMS=20097,16580 LABEL=d1

DISTANCE ATOMS=20096,16581 LABEL=d2

DISTANCE ATOMS=20111,2887 LABEL=d3

DISTANCE ATOMS=20535,2938 LABEL=d4

METAD ARG=d1     SIGMA=2    PACE=500  HEIGHT=1   FILE=HILLS LABEL=c1

ENERGY LABEL=ene

PRINT ARG=c1.bias STRIDE=500 FILE=COLVAR

PRINT ARG=ene FILE=COLVAR STRIDE=100

[Omar Valsson]

Hello, 

Your question is not clear, you are replying to an earlier message related to using ENERGY in LAMMPS while you say that you are using plumed2.2 in VMD. Can you clarify which MD code you are using? 

If you are using LAMMPS then you need to use a newer patch for LAMMPS that supports using the energy as CV, it is available in v2.3. The older patch available in v2.2 does not support the energy as CV. You will need to recompile LAMMPS after patching with the new patch. 

Regards,

Omar

--

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/a4540a7e-a5c3-4f87-bf8b-f1717a0d7bb3%40googlegroups.com.

[Omar Valsson]

Hello

Please keep the discussion on the mailing list for future reference. 

Ok, I understand that you are using the Plumed-GUI in VMD to analyze a trajectory generated with NAMD, correct? 

I do not have experience with Plumed-GUI in VMD but I cannot see how VMD can calculate the energy of the protein as it is not a MD code, therefore you cannot use the ENERGY CV as the energy is not available. You need to use that within the MD code (I am however not sure if the NAMD support the ENERGY CV). I am not too familiar with NAMD, but perhaps you can do a rerun on the trajectory. Furthermore, you should be aware that the ENERGY action will give the whole potential energy of the system which includes the solvent, is that what you are interested in? 

Omar

On 15 September 2017 at 14:54:53, João Fernandes (joaoro...@gmail.com) wrote:

I reported vmd because it is a vmd plugin, but in the case, a molecular dynamic generated by namd is used.

or I'm wrong, and it's not possible to do this program?

Enviado com Mailtrack

2017-09-15 9:49 GMT-03:00 Omar Valsson <omar.v...@gmail.com>:

Hello,

Ok, but I really don’t understand what you are trying to do. VMD is a molecular viewer and I do not understand how you can in there calculate the energy that needs to be passed to plumed. Do you perhaps mean NAMD? Please give further and clearer information, otherwise it is impossible to help you. 

Omar

On 15 September 2017 at 14:42:22, João Fernandes (joaoro...@gmail.com) wrote:

Sorry, it was my mistake.

my intention is to calculate the binding energy between 2 proteins.

I am using plumed2.2 in VMD.

I had the error of

so I made the correction in the file suggested here, so a new error came up.

to pass the total energy to plumed you should just add am instruction in the fix_plumed.cpp file.

After line:

p->cmd("setVirial",&virial[0][0]);

add:

p->cmd("setEnergy",&energy);

where energy is a double precision variable storing the energy.​

this is the new error

​terminate called after throwing an instance of 'PLMD::Exception'

  what():  

+++ Internal PLUMED error

+++ file Atoms.h, line 197

+++ message: assertion failed collectEnergy && energyHasBeenSet

 

I do not know if this modification also works for my problem, but it was the only solution I found for the initial error​

Enviado com Mailtrack

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/a4540a7e-a5c3-4f87-bf8b-f1717a0d7bb3%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

Att.
João Ronaldo Clemente Fernandes

Doutorando em Bioinformática - UFMG
Mestre em Biotecnologia Vegetal - UFLA
Bacharel em Biotecnologia - UNIPAC
skype jrcfernandes
e-mail joaoro...@gmail.com

linkedin http://br.linkedin.com/pub/joao-ronaldo-clemente-fernandes/4/396/773/
CV lattes http://lattes.cnpq.br/1677605981553830
+55 31 99216-3102
​+55 31 3223-2636​

--

Att.
João Ronaldo Clemente Fernandes

Doutorando em Bioinformática - UFMG
Mestre em Biotecnologia Vegetal - UFLA
Bacharel em Biotecnologia - UNIPAC
skype jrcfernandes
e-mail joaoro...@gmail.com

linkedin http://br.linkedin.com/pub/joao-ronaldo-clemente-fernandes/4/396/773/
CV lattes http://lattes.cnpq.br/1677605981553830
+55 31 99216-3102
​+55 31 3223-2636​

[João Fernandes]

I have a docking protein-protein result, and I was advised to use PLUMED-GUI to calculate the binding/interaction energy via metadynamics.

Since I already have a dynamic of this structure, I figured I could use the trajectory file (dcd) to get the results.

Enviado com Mailtrack

[Omar Valsson]

Hello, 

It does not make sense to use metadynamics on a already generated unbiased MD trajectory from NAMD, which I assume is what you are trying, as there is a feedback into the MD dynamics through the biasing forces that metadynamics generates. You will need to do a new NAMD run using PLUMED to do the metadynamics.

Omar

[João Fernandes]

Now I understand, I thought I could use an existing system dynamics for metadynamic.

Can you tell me a tutorial or script to perform this metadynamic?

[Giovanni Bussi]

Hi,

if you start from scratch, this is a good starting point:

https://en.wikipedia.org/wiki/Metadynamics

(it's the first match googling for "metadynamics").

If you navigate plumed tutorials you will find a lot of information

Giovanni

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