How to include vectors as CV in plumed

[Kartheekpitta(PhD), IIIT-H.]

Hi plumed users and experts,

I am trying to calculate free energy using well-tempered metadynamics for a protein containing two sub-units, choosing angle as collective variable (CV). The angle between two helices  which are present in two different subunits along  principal axis of each helices are considered as  CV . My question is as follows,

a) Can we consider CV as described above?

b) If yes, how to mention these vectors as CV in the plumed input file?

Thanks for a reply in advance.

[Liangzhen Zheng]

Hi Kartheekpitta,

I have the same problem. I am also thinking about using a helix as a vector. However, I didn't find much reference papers or materials about this.

And this kind of CV on a helix might distorted the Helix 2nd structure? It is another problem.

 

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

[Gareth Tribello]

Hello

I think this is a very interesting problem and although I don’t necessarily know how to solve it I do have a few suggestions.

Firstly, if you are able to define vectors that describe the orientation of your two helices you can certainly use the angle between these two vectors as a CV.  The simplest way I can think of to do this is to define a centres of mass based on the atoms at the two ends of the helix.  So if you have a helix involving atoms 1 to 100 say and a second helix involving atoms 101-200 you might proceed with an input something like this:

h1a: COM ATOMS=1-50 

h1b: COM ATOMS=51-100

h2a: COM ATOMS=101-150

h2b: COM ATOMS=150-200

a: ANGLE ATOMS=h1a,h1b,h2a,h2b

The idea here is that you are using the angle between the vectors connecting the centre of mass of the first 50 atoms in your helix to the centre of mass of the second 50 atoms in the helix.  I think that it is worth looking in the literature here though as there may well be better ways to define the orientation vector of an alpha helix.  The above idea is pretty crude.  

Secondly if you want to prevent the alpha helix from becoming distorted you might think about constraining the ALPHARMSD cv using UWALL.

Thanks

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 http://groups.google.com/group/plumed-users.
To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/a48f17d0-d77a-4728-a7e7-184e52a42969%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[joan.c...@gmail.com]

Hi everyone,

Just came across this thread and I am wondering if anything of the sort has been added since.

For context, I have a metric which is a vector with an arbitrary number of components, each of them ranging from 0 to 1. I am using the average of it as a CV and biasing it at 1 with a harmonic restraint, however, I have realised that even when 1 of this components is close to 1, if the rest of them aren't that component will still cause a strong bias because of the (x-x0) term in the derivative.  I have tried to define each component as an independent CV and bias them separately, but that make my calculations really slow.

Thanks

Joan

[Gareth Tribello]

HI 

If you use the latest version of PLUMED you can apply a bias on a vector of CVs.  If you look at the articles here:

It will give you some idea.  The short version is that, suppose I wasn to bias a bunch of distances.  I can do this like this:

d: DISTANCE ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6 

r: RESTRAINT ARG=d AT=0.1 KAPPA=0.3 

This input applies restraint forces on the distacnes between atoms 1 and 2, 3 and 4 and 5 and 6.  In other words, this input is equivalent to:

d1: DISTANCE ATOMS=1,2

r1: RESTRAINT ARG=d1 AT=0.1 KAPPA=0.3

d2: DISTANCE ATOMS=3,4

r2: RESTRAINT ARG=d2 AT=0.1 KAPPA=0.3

d3: DISTANCE ATOMS=5,6

r3: RESTRAINT ARG=d3 AT=0.1 KAPPA=0.3

However, the advantage of doing it the first way is that the calculation will be paralized over the distances.

I hope this helps

Gareth

To view this discussion visit https://groups.google.com/d/msgid/plumed-users/9142a37a-b2af-41d7-89f7-134fc6d369c9n%40googlegroups.com.

[joan.c...@gmail.com]

Hi Gareth,

Thanks for this, that's exactly what I am looking for.  Can this be used with any CV (including custom ones), and with CVs that take additional parameters other than ATOMS?

Context:

My CV takes a whole bunch of parameters:

GHOSTPROBE ...
LABEL=Psi0
PROBE_LABEL=0
ATOMS=Protein-H
PROBESTRIDE=25000 PRINT_PROTEIN
rmin=0.3076 DELTARMIN=0.0564
RMAX=0.4430 DELTARMAX=0.30
CMIN=0 DELTAC=1
PMIN=11.92 DELTAP=5.09
KPERT=1e-6 KXPLOR=1e-3 PERTSTRIDE=0
... GHOSTPROBE

So if I want to use 2 probes, it would become something like this:

Psi2: GHOSTPROBE PROBE_LABEL1=1 ATOMS1=Protein-H PROBESTRIDE1=25000 RMIN1=0.3076 DELTARMIN1=0.0564 RMAX1=0.4430 DELTARMAX1=0.30 CMIN1=0 DELTAC1=1 PMIN1=11.92 DELTAP1=5.09 KPERT1=1e-6 KXPLOR1=1e-3 PROBE_LABEL2=2 ATOMS2=Protein-H PROBESTRIDE2=25000 RMIN2=0.3076 DELTARMIN2=0.0564 RMAX2=0.4430 DELTARMAX2=0.30 CMIN2=0 DELTAC2=1 PMIN2=11.92 DELTAP2=5.09 KPERT2=1e-6 KXPLOR2=1e-3

But plumed (2.11) doesn't seem to understand the Psi2 line (works with your distance example though). I guess I could hardcode default parameters or read them from a text file, but I still need to specify PROBE_LABEL at runtime so that output files don't get all mixed up or overwritten.

Best,

Joan

[Gareth Tribello]

Hi Joan

You would need to do a little bit of coding to make the second input work.  I assume that your GHOSTPROBE action inherits from PLMD::Colvar?

If that is the case, then look at the code in src/colvar/Distance.cpp.  You should see that this has the following list of functions defined:

static void parseAtomList( const int& num, std::vector<AtomNumber>& t, ActionAtomistic* aa );

static unsigned getModeAndSetupValues( ActionWithValue* av );

static void calculateCV( const unsigned& mode, const std::vector<double>& masses, const std::vector<double>& charges,

                           const std::vector<Vector>& pos, std::vector<double>& vals, std::vector<std::vector<Vector> >& derivs,

                           std::vector<Tensor>& virial, const ActionAtomistic* aa );

The first two of these do things that were in the constructor for the old version of the command.  The last thing does a lot of what was done in the old version of calculate.

You will then notice that the command for registering the action is as follows:

typedef ColvarShortcut<Distance> DistanceShortcut;

PLUMED_REGISTER_ACTION(DistanceShortcut,"DISTANCE")

PLUMED_REGISTER_ACTION(Distance,"DISTANCE_SCALAR")

typedef MultiColvarTemplate<Distance> DistanceMulti;

PLUMED_REGISTER_ACTION(DistanceMulti,"DISTANCE_VECTOR")

The MultiColvarTemplate template class basically calls the calculateCv function function multiple times to calculate your vector of CVs.  Then once you have that vector of CVs you can pass it into the other actions.  It also makes use of the parseAtomList and getModeAndSetupValues to do the reading in of the atoms involved and working out what exactly the user wants to calculate (for the distance example this deals with the COMPONENTS and SCALED_COMPONENTS flags.  You may be able to do what is done in ANGLE for your version of the CV as there may not be options to do different things.)

Alternatively, if you don't want to change your code you can use two or more GHOSTPROBE commands and concatenate them together as follows.

cv1: GHOSTPROBE ...

cv2: GHOSTPROBE ...

vec: CONCATENATE ARG=cv1,cv2

 

You can then pass the object called vec to the restraint action.

I hope this helps

Gareth

To view this discussion visit https://groups.google.com/d/msgid/plumed-users/88af16af-199c-4649-96cd-c94c91323855n%40googlegroups.com.

[joan.c...@gmail.com]

Hi Gareth,

Thanks a lot for your help! Just before I modify the code though, does the CONCATENATE action make plumed calculate the CVs in parallel? (because if it does, that's exactly what I need)

Thanks

[Gareth Tribello]

No but if you do the coding the CVs will be computed in parallel.

Gareth

To view this discussion visit https://groups.google.com/d/msgid/plumed-users/9478a00f-057a-4c6a-9e26-d95d6da991f0n%40googlegroups.com.