Calculating Lechner & Dellago order parameters

[seanma...@gmail.com]

Hello, PLUMED team,

I'm interested in using PLUMED's crystallization module to implement variants of like the bond orientational order parameters (OPs) devised by Lechner and Dellago, which I will denote here by q_LD(i) for the OP associated with the i-th atom. For those not familiar with it already, q_LD(i) is constructed by first taking the average of the Steinhardt vectors q_l,m(i) of the neighbors of i and i itself, and then taking the complex norm. In water, for example, this means averaging the q_l,m(i) of 5 atoms on average (4 neighbors plus the "central" water itself).

I've been studying the code for a few days now, and I'm impressed by its versatility. In that spirit, I'd like to take advantage of the features already contained in PLUMED. I have a few questions about the inner workings of the code that I hope you can answer.

1. By carefully choosing D_MAX in the switching function, can I guarantee that any atom "i" will have access to the positions to all of its neighbors in the first and second hydration shell? My impression is that this is how the LinkCells parallelizes the calculation of quantities like q_l,m(i).
   
   
2. For each processor in a simulation, does PLUMED (for each rank) only grab the atoms local to that processor? Or is there a global atom collection that takes place?
   
1. I'm concerned about performance because I'd eventually like to calculate q_LD(i) for fairly large systems (hundreds to thousands of atoms)
   
   
3. Is it possible to calculate the average value of a vector of multicolvars in a probe volume of arbitrary geometry (e.g. spheres, cylinders, etc.)?
   
1. Is it possible to then bias this average value?

It it's useful to know, my ultimate goal is to study crystallization in water using GROMACS.

Sincerely,

Sean Marks

Patel Group

Penn CBE

[Gareth Tribello]

Hello

I don’t know if you are already aware but we have just published a paper on the crystallisation stuff that we have been implementing in the crystallisation module.  It is available here:

http://pubs.acs.org/doi/abs/10.1021/acs.jctc.6b01073

This might be useful in terms of getting up and running.  In terms of your questions.  

1. I think the answer to this is yes.  For any value of d>d_max the switching function is zero and thus we do not need to calculate it.  PLUMED thus uses D_MAX as the link cells cutoff.  Incidentally, (and I’m sorry if I am stating the obvious) you should always use a D_MAX parameter as otherwise link cells are not used and the code will thus be impossibly slow.

2. Is this a question about domain decomposition?  PLUMED gathers the positions of all the atoms it needs from the underlying MD engine and then passes them to all nodes.  Everything within PLUMED is thus done with particle decomposition.  My understanding is that the reason for this is transferability.  Different MD codes having different ways to doing the domain decomposition algorithm so it is difficult to work out a general way to PLUMED work with all domain decomposed codes.  On the plus side though if you look at that paper we were looking at some pretty massive systems.  Having said that though we were only really doing post processing.

3. Yes you do something like this:

q6: Q6 SPECIES=1-100 SWITCH={RATIONAL D_0=1 R_0=1 D_MAX=4}

sp: INSPHERE ATOM=101 SPECIES=q6 RADIUS={RATIONAL D_0=10 R_0=3 D_MAX=15}  MEAN

PRINT ARG=sp.mean FILE=colvar

This is calculating the average Q6 value inside a sphere of radius 10 nm centered on the position of atom 101.  Notice there are some other volumes to find out more you can look at the values of CVs in as detailed here:

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

You can also implement new volumes relatively easy.

To do the Lechner and Dellago order parameters you need something like this:

q6: Q6 SPECIES=1-100 SWITCH={RATIONAL D_0=1 R_0=1 D_MAX=4}

q6l: LOCAL_AVERAGE SPECIES=q6 SWITCH={RATIONAL D_0=1 R_0=1 D_MAX=4}

This is doing a local average (of the complex vectors) in a sphere around each of the atoms input to q6.  You should also be able to feed the output from q6l into INSPHERE so that you can calculate the values of these quantities in a shape.

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/72fc39f3-5a4d-4ad6-b8d6-2891b61c009c%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Sean Marks]

Hello, Gareth,

Thank you for such a prompt and helpful reply! I've read the PLUMED 2 paper, but I didn't realize there was another on crystallization; I will make sure to read it thoroughly ASAP. I have just one follow-up question: Once I have the Lechner & Dellago parameters for the atoms in the sphere, is it possible to bias their average value, or some other function of the atom-wise quantities?

Thanks for your time!

Sean

[Gareth Tribello]

Hello

No problem.  It should be possible to bias the average value of the Dellago parameters or anything you can calculate with multicolvar.  I think you would be the first to do it though.  When you are doing something newish like this it is worth having a check of the derivatives before running anything.  I did test them and they should be right but it does no harm to double check,

You can do this check using driver.   Suppose that you have a PLUMED input with a restraint i.e. an input that is definitely going to apply a force.  An input like this would be an example:

d1: DISTANCE ATOMS=1,2

RESTRAINT ARG=d1 KAPPA=10 AT=1.0

The important thing is that there need to be a restraint in there otherwise no force is applied.  

If you then run driver on a short trajectory using something like the following command:

plumed driver —ixyz trajectory.xyz —debug-forces forces.num

Analytical and numerical derivatives of the forces will be output in a file called forces.num.  If you could let me know if there is a problem I would greatly appreciate it.

Have fun

Gareth

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/4b9fd575-e60e-4248-b925-96da47a117a2%40googlegroups.com.

[Sean Marks]

Will do. Thanks again for all your help!

[Sean Marks]

Hi, Gareth,

How might I use/upgrade PLUMED to calculate a linear combination of Lechner and Dellago's order parameters in a particular volume? I'm thinking about a new data member for INSPHERE/INCYLINDER like "MEAN", which can then be acted upon by other PLUMED actions (e.g. analyzed, biased).

From studying the code, it seems that the keyword MEAN is implemented as a FunctionVessel. Would my desired result be as simple as implementing a new FunctionVessel in the spirit of MEAN, with keywords registered with INSPHERE? Or are there more complex dependencies (like ActionWithVessel, VesselRegister, etc.) that need to be considered?

As an aside: What is a "vessel," exactly? Is it an interface which makes it relatively easy to parallelize multicolvar calculations, with a large buffer of doubles that can be passed around easily by MPI?

Best,

Sean

[Gareth Tribello]

Hello Sean

Thanks for your email.  In terms of your question it would depend on what linear combination you wanted to calculate.  It might be possible to do it by implementing a new class that inherited from FunctionVessel.  But I would need to know something more about the coefficients in order to say something more certain.  

In terms of what a vessel is your definition is pretty good.  Essentially in multicolvar you accumulate quantities in a long std::vector.  Then at the end of the parallel loop you do an all gather on this std::vector.  The vessels give you an easy interface that allows you to extract what you need from this std::vector in terms of values and derivatives for the final quantities you are trying to compute. 

Thanks

Gareth

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/23021158-e238-404c-b86d-2ac9582dd1d9%40googlegroups.com.

[Sean Marks]

Thanks for clarifying the inner workings of PLUMED for me.

To be more precise, what I would like to do is calculate the sum of the Lechner and Dellago order parameters, and then divide that sum by a constant (which could be hard-coded). I then want to bias the resulting sum. Though the constant is unimportant, because it could be factored out and included as part of the definition of the biasing parameters, KAPPA and SLOPE.

[Gareth Tribello]

I think you can do this without modifying the code.  Just do:

q6: Q6 SPECIES=1-100 …

q6a: LOCAL_AVERAGE SPECIES=q6 … 

vv: INSPHERE DATA=q6a SUM ...

COMBINE ARG=vv.sum COEFFICIENTS=<constant>

My constant here would be one over your constant as the coefficient in combine is a multiplicative constant.

The sum option might not be in INSPHERE but you can add it easily by justing put:

keys.use(“SUM”);

In the registerKeywords function for the action where you want this.  

Gareth

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/6a9fdd92-107a-456c-b57e-14ad4d51410f%40googlegroups.com.

[Sean Marks]

Hi, Gareth,

That's fantastic!

Right now, I'm incrementally testing my own force calculations by computing forces for simpler order parameters (L&D's OP gives rise to some pretty nasty derivatives!). I'm currently looking to test PLUMED vs. my own code by computing the classic probe-volume-average Steinhardt parameter Q_6 (Q_l, below)

Does INSPHERE's VMEAN option compute VMEAN(Q6) as I've defined it above when you pass it a Q6 multicolvar using the DATA keyword? Or does it compute some other variant?

Thanks!

Sean

PS: I've attached a PNG of the equations in case the image I tried to include above doesn't render properly.

[Sean Marks]

One more thing: When I used the --debug-forces option on a single frame in a trajectory with 16,500 atoms, the resulting file had 49,509 = 3*16,500 + 9 entries. Are these last 9 entries the components of the virial? Here's my input file:

# Water oxygens (TIP4P/Ice)
ow: GROUP ATOMS=1-16500:4

q6: Q6 SPECIES=ow SWITCH={GAUSSIAN R_0=0.02 D_0=0.32 D_MAX=0.34}

center: FIXEDATOM AT=2.5,2.5,2.5

sphere: INSPHERE ATOM=center DATA=q6 RADIUS={GAUSSIAN D_0=0.50 R_0=0.02 D_MAX=0.52} VMEAN

restraint: RESTRAINT ARG=sphere.vmean AT=0.5 KAPPA=1000.0 SLOPE=0.0

PRINT ARG=sphere.vmean,restraint.bias,restraint.force2 FILE=bias_Q6.out

Also, the numerical derivatives column seems to contain garbage values (massive numbers). Does --debug-forces only compute the analytic derivatives by default?

[Gareth Tribello]

Hello

As a first comment please define your symbols when you write equations.  It is very hard to answer if you don’t explain the symbols.

I think, however, that this is not what VMEAN is computing if you pass to INSPHERE.  Lets consider the following input:

q6: Q6 SPECIES=1-100 SWITCH={RATIONAL R_0=0.1} VMEAN

PRINT ARG=q6.* FILE=colvar

Using your nomenclature this computes something like:

[Sean Marks]

Hi, Gareth,

First of all, I apologize: You're absolutely correct, I should thoroughly define all of my symbols. I've been staring at them for so long that I made the mistake of assuming what was obvious to me would be obvious to someone else looking at them.

Attached is a PDF with what I hope are better-defined symbols. I attempted to use the nomenclature of PLUMED's manual wherever possible, and I added primes to the new quantities I'm defining. Could you confirm whether I understand what PLUMED is doing?

Thanks for your time!

Sean

[Sean Marks]

And here is the PDF I meant to include!

[Gareth Tribello]

Hello

In response to your question. If you use:

q6: Q6 SPECIES=1-100 SWITCH={RATIONAL R_0=1.0} MEAN

It calculates the average of the norm of the vector.

If you use:

q6: Q6 SPECIES=1-100 SWITCH={RATIONAL R_0=1.0} VMEAN

It calculates the average vector and then calculates the norm of this average vector.

Everything else you do, be putting it through INSPHERE, using LOCAL_AVERAGE, etc, just replaces the above averages with weighted averages. In these different quantities the weights are calculated in different ways.

I think from your maths that what you want to bias is the quantity computed by Q6 + INSPHERE + VMEAN. The only thing is that PLUMED is missing a factor of \sqrt{ 4*\pi / (2l+1) }. I don’t think this matters though as it is just a multiplicative constant.

I hope this helps
Gareth

> On 23 Feb 2017, at 21:33, Sean Marks <seanma...@gmail.com> wrote:
>
> <PLUMED_VMEAN.pdf>

[Sean Marks]

Hi, Gareth,

Thanks for your help! I decided to first test whether my code agreed with PLUMED for calculating the value of Q6-INSPHERE-VMEAN, before going through and checking forces. Unfortunately, I'm getting some "nan" values for the order parameter when I use 'driver' to post-process a simulation of liquid TIP4P/Ice in a box I performed using unmodified GROMACS. Here is my plumed.dat:

# Water oxygens: atoms 1-N with stride 4 (TIP4P/Ice)
ow: GROUP ATOMS=1-16500:4

# Calculates q_{l,m}(i)
q6: Q6 SPECIES=ow SWITCH={GAUSSIAN D_0=0.32 R_0=0.01 D_MAX=0.34}

# Fixed virtual atom which serves as the probe volume's center (pos. in nm)
center: FIXEDATOM AT=2.5,2.5,2.5

# Probe volume
sphere: INSPHERE ATOM=center DATA=q6 RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52} VMEAN

# Include factor of sqrt(4*pi/(2*l+1) (for l=6) to get Q_6,v (Steinhardt definition)
q6v: COMBINE ARG=sphere.vmean COEFFICIENTS=0.98318 PARAMETERS=0.0 POWERS=1.0 PERIODIC=NO

# Bias the mean value of the OP in the sphere
restraint: RESTRAINT ARG=q6v AT=0.5 KAPPA=1000.0 SLOPE=0.0

print: PRINT ARG=q6v,restraint.bias FILE=bias_Q6_plumed.out

And here is the command I use to run driver:

mpi-plumed_2.3.0 driver --plumed 'plumed.dat' --timestep 0.002 --trajectory-stride 500 --mf_xtc 'myxtcfile.xtc'

I only see 'nan' for 3 out of about 10,000 frames; the rest of the CV values from PLUMED agree with my own calculations to within <1%. My own code doesn't calculate any 'nans'. Do you see anything wrong with my input? Or could this perhaps be a bug in PLUMED? I'm happy to provide more information.

Best,
Sean

[Gareth Tribello]

Hi Sean

My guess would be that there are no atoms in the volume of interest and hence the nan.  Can you send me a couple of frames that give nans so I can take a look.  If it is what I think it should be easy to fix.

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 https://groups.google.com/group/plumed-users.

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/91626f94-5ff1-498b-9cc7-b78ca9cb9e92%40googlegroups.com.

[Sean Marks]

Hi, Gareth,

Here are the three frames that given nan. In each case, my code says there are at least 10 applicable atoms (i.e. water oxygens) in the volume, so I don't think it's an issue with the frames themselves. It bears mentioning that with my output, PLUMED prints some potentially incorrect information to stdout:

PLUMED: Action Q6
PLUMED:   with label q6
PLUMED:   Steinhardt parameter of central atom and those within 0.01.  Using gaussian swiching function with parameters d0=0.32
PLUMED:   keyword SPECIES takes atoms : 1 5 9 ...
PLUMED: Action FIXEDATOM
PLUMED:   with label center
PLUMED:   serial associated to this virtual atom is 16501
PLUMED:   AT position 2.5 2.5 2.5
PLUMED: Action INSPHERE
PLUMED:   with label sphere
PLUMED:   calculating q6 inside region of insterest
PLUMED:   added component to this action:  sphere.vmean
PLUMED:   value sphere.vmean contains the norm of the mean vector
PLUMED:   center of sphere is at position of atom : 16501
PLUMED:   radius of sphere is given by 0.01.  Using gaussian swiching function with parameters d0=0.5
PLUMED: Action COMBINE
PLUMED:   with label q6v
PLUMED:   with arguments sphere.vmean
PLUMED:   with coefficients: 0.983180
PLUMED:   with parameters: 0.000000
PLUMED:   and powers: 1.000000
PLUMED: Action RESTRAINT
PLUMED:   with label restraint
PLUMED:   with arguments q6v
PLUMED:   added component to this action:  restraint.bias
PLUMED:   at 0.500000
PLUMED:   with harmonic force constant 1000.000000
PLUMED:   and linear force constant 0.000000
PLUMED:   added component to this action:  restraint.force2
PLUMED: Action PRINT
PLUMED:   with label print
PLUMED:   with stride 1
PLUMED:   with arguments q6v restraint.bias
PLUMED:   on file bias_Q6_plumed.out
PLUMED:   with format  %f

Under the Q6 action, I'm not sure what is meant by "and those within 0.01"; is this just a way to say "r_0 in the switching function is 0.01 nm"? Or is it actually doing some kind of local average? My goal is to have a neighbor sphere with a radius of 0.32 nm, with a switching function that goes to zero at d_max = 0.34. I have the same question about what I've highlighted in the INSPHERE action, where my goal is a probe sphere of radius 0.5 nm.

On a random note, in the online documentation for INSPHERE in v2.3.0, the definition for w(x_i, y_i, z_i) seems to be missing. And in the line below that, I believe an "f" is missing from "f(s_i)".

Best,
Sean

[Gareth Tribello]

Hello again Sean

Thanks for these frames.  I had a run of them and I don’t get nan, which is weird.  I am running with the latest master version of plumed which you can download using:

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

Can you have a go with that version on your machine to see if it is something machine specific or if it is just that there is an error in an old version of plumed that has now been repaired.

Separately on the switching function it is just that the output command reads:

std::cout<<“and those within “<<r_0<<“.  Using gaussian switching function with parameters d0=“<<0.23<<std::endl;

which I admit is perhaps misleading for a log message.

Thanks

Gareth

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/66c72c35-1eab-4d23-9508-fab050ceead9%40googlegroups.com.

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

<traj_frame_1177ps.xtc><traj_frame_2951ps.xtc><traj_frame_7264ps.xtc>

[Sean Marks]

Hi, Gareth,

Thanks for working with me to resolve this issue. I installed the 2.4.0-dev version using the git repository you indicated and reran the calculations, but I'm still getting NaN on the same frames. What do you think might be the issue? Perhaps I'm compiling/linking incorrectly; here's the command I ran to configure PLUMED:

./configure --program-suffix=_git_mpi \
            --enable-modules=crystallization \
            LDFLAGS="-L/usr/local/lib" \
            CPPFLAGS="-I/usr/local/include" \
            LIBS="-lxdrfile"

I've also attached my config.log. I haven't done the regtests for 2.4.0-dev, but my 2.3.0 installation passed all applicable tests.

Sean

[Sean Marks]

By the way, I should mention I'm using an iMac Retina 5k running El Capitan (v10.11.16) with an Intel Core i5.

[Gareth Tribello]

Hello again

Are you running with MPI/Open_mp when you do the calculations or are you just running on one node?  I ran the calculations on one node so maybe it breaks if it runs in parallel?

Gareth

On 08 Mar 2017, at 21:56, Sean Marks <seanma...@gmail.com> wrote:

Hi, Gareth,

Thanks for working with me to resolve this issue. I installed the 2.4.0-dev version using the git repository you indicated and reran the calculations, but I'm still getting NaN on the same frames. What do you think might be the issue? Perhaps I'm compiling/linking incorrectly; here's the command I ran to configure PLUMED:

./configure --program-suffix=_git_mpi \
            --enable-modules=crystallization \
            LDFLAGS="-L/usr/local/lib" \
            CPPFLAGS="-I/usr/local/include" \
            LIBS="-lxdrfile"


I've also attached my config.log. I haven't done the regtests for 2.4.0-dev, but my 2.3.0 installation passed all applicable tests.

Sean

On Wednesday, March 8, 2017 at 4:05:32 PM UTC-5, Gareth Tribello wrote:

Hello again Sean

Thanks for these frames.  I had a run of them and I don’t get nan, which is weird.  I am running with the latest master version of plumed which you can download using:

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

Can you have a go with that version on your machine to see if it is something machine specific or if it is just that there is an error in an old version of plumed that has now been repaired.

Separately on the switching function it is just that the output command reads:

std::cout<<“and those within “<<r_0<<“.  Using gaussian switching function with parameters d0=“<<0.23<<std::endl;

which I admit is perhaps misleading for a log message.

Thanks

Gareth

On 07 Mar 2017, at 17:23, Sean Marks <seanma...@gmail.com> wrote:

Hi, Gareth,

Here are the three frames that given nan. In each case, my code says there are at least 10 applicable atoms (i.e. water oxygens) in the volume, so I don't think it's an issue with the frames themselves. It bears mentioning that with my output, PLUMED prints some potentially incorrect information to stdout:

PLUMED: Action Q6
PLUMED:   with label q6


PLUMED:   Steinhardt parameter of central atom and those within 0.01.  Using gaussian swiching function withparameters d0=0.32

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/132d5280-1371-41f5-b3da-74439828c666%40googlegroups.com.

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

<config.log>

[Sean Marks]

I believe I've been running it on only one processor, though I did compile it with MPI enabled. The messages printed to stdout say that I'm only using 1 thread on 1 node.

When I do run it with mpirun on 4 cores, I get NaNs for the same frames.

[Sean Marks]

Now this is odd. Up until now, I've been running PLUMED and my own code on the xtc file of the whole trajectory, and I see NaNs. But when I run it on the individual frames I sent you, I get the following finite answers:

time [ps]      q6v               restraint.bias
1177            0.154901      59.546667
2951            0.160024      57.791988
7264            0.100374      79.850323

To view this discussion on the web visit <a href="https://groups.google.com/d/msgid/plumed-users/132d5280-1371-41f5-b3da-74439828c666%40googlegroups

[Sean Marks]

Would you be able to try running the same calculation on the attached set of 20 frames? They're in the vicinity of 1177 ps (range 1170-1190 ps). If I run t = 1177 ps on its own, PLUMED gives me a finite answer (albeit one which is very different from what my code computes). But if I run PLUMED with either the 20-frame sub-trajectory or the whole trajectory, I get NaN each time.

For reference, here is what my code computes for these three frames. Generally, my results are the same as PLUMED to within 1% (at most 2.5%).

time [ps]      q6v               restraint.bias       q6v(my code)
1177            0.154901      59.546667          0.135736
2951            0.160024      57.791988          0.107751
7264            0.100374      79.850323          0.19153

[Gareth Tribello]

Hello again

I ran the trajectory you sent with PLUMED and the following command:

plumed driver --mf_xtc near_1177ps.xtc

The input I used was this one:

 Water oxygens: atoms 1-N with stride 4 (TIP4P/Ice)

ow: GROUP ATOMS=1-16500:4

# Calculates q_{l,m}(i)

q6: Q6 SPECIES=ow SWITCH={GAUSSIAN D_0=0.32 R_0=0.01 D_MAX=0.34}

# Fixed virtual atom which serves as the probe volume's center (pos. in nm)

center: FIXEDATOM AT=2.5,2.5,2.5

# Probe volume

sphere: INSPHERE ATOM=center DATA=q6 RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52} VMEAN

# Include factor of sqrt(4*pi/(2*l+1) (for l=6) to get Q_6,v (Steinhardt definition)

q6v: COMBINE ARG=sphere.vmean COEFFICIENTS=0.98318 PARAMETERS=0.0 POWERS=1.0 PERIODIC=NO

# Bias the mean value of the OP in the sphere

restraint: RESTRAINT ARG=q6v AT=0.5 KAPPA=1000.0 SLOPE=0.0

print: PRINT ARG=q6v,restraint.bias FILE=bias_Q6_plumed.out

And the output was this:

! FIELDS time q6v restraint.bias

 0.000000 0.144393 63.228105

 1.000000 0.127395 69.417253

 2.000000 0.093932 82.445596

 3.000000 0.147839 62.008659

 4.000000 0.144977 63.020582

 5.000000 0.144423 63.217472

 6.000000 0.092486 83.033848

 7.000000 0.135736 66.344089

 8.000000 0.174318 53.034359

 9.000000 0.247137 31.969771

 10.000000 0.141759 64.168363

 11.000000 0.099090 80.364310

 12.000000 0.127046 69.547440

 13.000000 0.172266 53.704776

 14.000000 0.147718 62.051421

 15.000000 0.144188 63.301036

 16.000000 0.149100 61.565573

 17.000000 0.161993 57.124459

 18.000000 0.097169 81.136386

 19.000000 0.155983 59.173852

 20.000000 0.104583 78.177207

As you can see I didn’t get any nan's.  Could there be something different in your version of molfile? This is what PLUMED is using to read in the trajectory.

In terms of the three frames you sent yesterday.  I get exactly the same numbers you do when I run PLUMED. 

Curiouser and curiouser.

Gareth

On 08 Mar 2017, at 23:08, Sean Marks <seanma...@gmail.com> wrote:

<near_1177ps.xtc>

[Giovanni Bussi]

Hi guys,

if there is a problem with xtc reader, please compare

--mf_xtc

with

--ixtc (this needs xdrfile library to be compiled).

I tried both options with the following plumed.dat file:

DUMPATOMS ATOMS=1-16500 FILE=test-mf.gro

and I get two slightly different box shapes:

16503c16503

<    5.0091825    5.0091825    5.0091825    0.0000000    0.0000000    0.0000000    0.0000000    0.0000000    0.0000000

---

>    5.0091824    5.0091824    5.0091822    0.0000000    0.0000000    0.0000000    0.0000000    0.0000000    0.0000000

33006c33006

<    5.0157738    5.0157738    5.0157738    0.0000000    0.0000000    0.0000000    0.0000000    0.0000000    0.0000000

---

>    5.0157738    5.0157738    5.0157737    0.0000000    0.0000000    0.0000000    0.0000000    0.0000000    0.0000000

Notice that box read with mf plugin (second) is not exactly cubic, whereas box read with ixtc library (first) is exactly cubic.

Gareth: perhaps there is a bug in the code that only appears when the box is exactly cubic?

Giovanni

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/7cca06f7-7609-400e-ae9c-2b23b4124ac6%40googlegroups.com.

[Sean Marks]

Hi, team,

Gareth, I get exactly the same numbers as you do except for t = 7 ps (originally 1177 ps) when I use mf_xtc. When I use ixtc with the same installation of the xdrfile library I link to my own program, I get:

#! FIELDS time q6v restraint.bias
 0.000000 0.144393 63.228119
 1.000000 0.127395 69.417232
 2.000000 0.093932 82.445490
 3.000000 0.147839 62.008660
 4.000000 0.144977 63.020641
 5.000000 0.144423 63.217521
 6.000000 0.092486 83.033794
 7.000000 nan nan
 8.000000 0.174318 53.034346
 9.000000 0.247137 31.969790
 10.000000 0.141759 64.168392
 11.000000 0.099090 80.364279
 12.000000 0.127046 69.547494
 13.000000 0.172266 53.704756
 14.000000 0.147717 62.051531
 15.000000 0.144188 63.301068
 16.000000 0.149100 61.565547
 17.000000 0.161993 57.124447
 18.000000 0.097169 81.136330
 19.000000 0.155983 59.173815
 20.000000 0.104583 78.177226

The values of "q6v" are exactly the same between my mf_xtc and ixtc calculations, and restraint.bias is the same up to the 3rd decimal place.

Sean

To view this discussion on the web visit <a href="<a href="https://groups.google.com/d/msgid/plum

[Sean Marks]

I used GROMACS' trjconv to convert the 20-frame xtc file to both a gro file and a PDB file (see https://drive.google.com/open?id=0Bw86zVjv5VTRMkkwVGhqMldRbXc). I ran again using --igro and --mf_pdb, and I'm getting the same behavior. I think this rules out problems with the xtc readers themselves.

Sean

<span style="font-family:Helvetica;font-size:12px;font-style:normal;font-var

[Gareth Tribello]

Hello again Sean

I downloaded your pdb and gro files from that google drive and ran PLUMED on them and I still don’t get any nans.  Were you getting nans when you ran PLUMED on those trajectories?  I haven’t tried with --ixtc yet

Thanks

Gareth

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/39365f5b-dd4d-470f-ade6-d5e997a53c22%40googlegroups.com.

[Sean Marks]

Hi, Gareth,

Yes, I am. I'm baffled as to why.

Sean

> On 23 Feb 2017, at 21:33, Sean Marks <seanma...@<a href="http://gmail.com/" rel="nofollow" target="_blank" onmousedown="this.href='http://gmail.com/';return true;" onclick="this.href='http://

[Sean Marks]

Hi again,

If you move the center of the probe volume from (2.5, 2.5, 2.5) by a few Angstroms, PLUMED produces finite values for me. This makes me think the issue might be somewhere in INSPHERE.

Sean

[Sean Marks]

Hi, PLUMED team,

I've been looking at the implementation of the Steinhardt class in the source code, and I think there may be several mistakes in the calculation of the spherical harmonics and their derivatives. For one, I believe the calculation of the associated Legendre polynomial P_{l,m}(x) leaves out the factor of (1 - x^2)^{m/2} (see http://mathworld.wolfram.com/AssociatedLegendrePolynomial.html). I also see potential problems with the derivatives of "z" = exp(i*m*phi) with respect to {x,y,z}: I think there are errors in the chain rule here.

I have worked out all of the calculations myself, and I've attached them for your reference. I can't guarantee that they're without error themselves, but I have checked my work numerous times and I'm fairly confident in the results. I could also type up my derivation of the derivatives of the harmonics, if that helps.

Best,

Sean

[Sean Marks]

PS: My write-up has a fairly pedantic tone because I wrote it as a tutorial for members of my research group.