Convert bonded potential to generic table

[Vitalie Botan]

Hi there,

I have a Lammps trajectory and corresponding topology and mapping XML files, so everything seems to work fine with csg_boltzmann and csg_stat scripts and I am able to generate CG-tables both for bonded and non-bonded interactions. However, when I try to convert  resampled and extrapolated bonded potential with the following script csg_call convert_potential lammps

got this message:

ERROR:                                                                                              

potential_to_generic.sh: conversion of bonded interaction to generic tables is not implemented yet!

Is there any proper way to prepare a smoothed extrapolated bonded potential for Lammps? I've tried csg_resample with --derivative option together with the smoothed potential from csg_boltzmann output, but all the table values were marked with "o" sign, meaning they will be extrapolated later, and exactly this scenario one has to avoid when preparing tabulated potential for Lammps.

--Vitalie

PS: I am using 1.3-dev version of votca.

[Christoph Junghans]

2014-12-01 18:11 GMT-07:00 Vitalie Botan <vitali...@gmail.com>:
> Hi there,
>
> I have a Lammps trajectory and corresponding topology and mapping XML files,
> so everything seems to work fine with csg_boltzmann and csg_stat scripts and
> I am able to generate CG-tables both for bonded and non-bonded interactions.
> However, when I try to convert resampled and extrapolated bonded potential
> with the following script csg_call convert_potential lammps
> got this message:
>
> ERROR:
>
> potential_to_generic.sh: conversion of bonded interaction to generic tables
> is not implemented yet!

Most of the code is already there, it is just not enabled, because I
have never used tabulated bonded interaction in LAMMPS.
Do you have a specification for LAMMPS' bonded tables? (i.e. what's
in the 1st, 2nd, 3rd column for bond, angle and dihedral tables?)

>
>
> Is there any proper way to prepare a smoothed extrapolated bonded potential
> for Lammps?

You could use csg_call convert_potential gromacs instead, gromacs'
bonded tables are just three columns. (r, U(r), F(r))

I've tried csg_resample with --derivative option together with
> the smoothed potential from csg_boltzmann output, but all the table values
> were marked with "o" sign, meaning they will be extrapolated later, and
> exactly this scenario one has to avoid when preparing tabulated potential
> for Lammps.

You could change it by hand from "o" to "i", or run the
csg_call potential extrapolate [OPTIONS] input output

Christoph

>
> --Vitalie
>
>
> PS: I am using 1.3-dev version of votca.
>

> --
> You received this message because you are subscribed to the Google Groups
> "votca" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to votca+un...@googlegroups.com.
> To post to this group, send email to vo...@googlegroups.com.
> Visit this group at http://groups.google.com/group/votca.
> For more options, visit https://groups.google.com/d/optout.



--
Christoph Junghans
Web: http://www.compphys.de

[Vitalie Botan]

On Tuesday, December 2, 2014 4:47:22 PM UTC+1, Christoph Junghans wrote:

2014-12-01 18:11 GMT-07:00 Vitalie Botan <vitali...@gmail.com>:

Most of the code is already there, it is just not enabled, because I
have never used tabulated bonded interaction in LAMMPS.
Do you have a specification for LAMMPS'  bonded tables? (i.e. what's
in the 1st, 2nd, 3rd column for bond, angle and dihedral tables?)

Hi Christoph,

Thanks for the prompt reply, it seems I overlooked the fact, that gromacs version of the script produces not only xvg-table, but a normal one as well.

Regarding the format of the bonded tabulated potential, it is very similar to the non-bonded format (which is already implemented in votca) with few exceptions I have copied from the Lammps manual below:

1) A section begins with a non-blank line whose 1st character is not a "#"; blank lines or lines starting with "#" can be used as comments between sections. The first line begins with a keyword which identifies the section. The next line lists (in any order) one or more parameters for the table. Each parameter is a keyword followed by one or more numeric values. The parameter "N" is required and its value is the number of table entries that follow. Following a blank line, the next N lines list the tabulated values. On each line, the 1st value is the index from 1 to N, the 2nd value is the bond length r (in distance units), the 3rd value is the energy (in energy units), and the 4th is the force (in force units). The bond lengths must range from a LO value to a HI value, and increase from one line to the next. If the actual bond length is ever smaller than the LO value or larger than the HI value, then the bond energy and force is evaluated as if the bond were the LO or HI length.

2) For angle potential, on each line, the 1st value is the index from 1 to N, the 2nd value is the angle value (in degrees), the 3rd value is the energy (in energy units), and the 4th is -dE/d(theta) (also in energy units). The 3rd term is the energy of the 3-atom configuration for the specified angle. The last term is the derivative of the energy with respect to the angle (in degrees, not radians). Thus the units of the last term are still energy, not force. The angle values must increase from one line to the next. The angle values must also begin with 0.0 and end with 180.0, i.e. span the full range of possible angles.

3) For dihedral potential, on each line, the 1st value is the index from 1 to N, the 2nd value is the angle value, the 3rd value is the energy (in energy units), and the 4th is -dE/d(phi) also in energy units). The 3rd term is the energy of the 4-atom configuration for the specified angle. The 4th term (when present) is the negative derivative of the energy with respect to the angle (in degrees, or radians depending on whether the user selected DEGREES or RADIANS). Thus the units of the last term are still energy, not force. The dihedral angle values must increase from one line to the next. 

Dihedral table splines are cyclic. There is no discontinuity at 180 degrees (or at any other angle). In general, the first angle in the list can have any value (positive, zero, or negative). However the range of angles represented in the table must be strictly less than 360 degrees (2pi radians) to avoid angle overlap. (You may not supply entries in the table for both 180 and -180, for example.) If the user's table covers only a narrow range of dihedral angles, strange numerical behavior can occur in the large remaining gap.

Hope these restrictions don't pose any problem for implementing in the script. Thanks in advance.

Best,

--Vitalie

[Christoph Junghans]

I hacked something together quickly:
<https://code.google.com/p/votca/source/detail?r=1e124eb53c82c3fd451e08f57065ac10ec17eb4f&repo=csg>
Dihedral tables aren't support in their full beauty yet, but also for
Gromacs the periodicity of the table is not imposed.

Christoph

>
> Best,
>
> --Vitalie

[Vitalie Botan]

Thanks a lot, Christoph. Works like a charm!

Vitalie

On Tuesday, December 2, 2014 11:17:10 PM UTC+1, Christoph Junghans wrote:

2014-12-02 9:53 GMT-07:00 Vitalie Botan <vitali...@gmail.com>:
I hacked something together quickly:
<https://code.google.com/p/votca/source/detail?r=1e124eb53c82c3fd451e08f57065ac10ec17eb4f&repo=csg>
Dihedral tables aren't support in their full beauty yet, but also for
Gromacs the periodicity of the table is not imposed.

Christoph

[Vitalie Botan]

Hi Christoph,

It looks there is a problem with the angle potential routine, it fails with the following error:

ERROR: 

do_external: subscript

/opt/votca/share/votca/scripts/inverse/table_extrapolate.pl --function linear --avgpoints 3 --region leftright ABB.pot.smooth.lQm2b ABB.pot.extrapol.CUPup

(from tags table extrapolate) failed

For details see the logfile

I searched the log file and found a command, which could cause this error: 

Running subscript 'potential_to_lammps.sh ABB.extra.pot ABB.bonded.pot' (from tags convert_potential lammps) dir /opt/votca/share/votca/scripts/inverse

Convert ABB.extra.pot to ABB.bonded.pot

Running critical command 'sed -n s/.*cg_bonded\.\([^[:space:]]*\) .*/\1/p'

Running critical command 'mktemp ABB.pot.scale.XXXXX'

Running subscript 'table_linearop.pl --on-x ABB.extra.pot ABB.pot.scale.d4RZ0 57.2957795 0' (from tags table linearop) dir /opt/votca/share/votca/scripts/inverse

table_linearop.pl: ABB.extra.pot to ABB.pot.scale.d4RZ0 with x' = 57.2957795*x + 0 

Running critical command 'mktemp ABB.pot.smooth.XXXXX'

Running critical command 'csg_resample --in ABB.pot.scale.d4RZ0 --out ABB.pot.smooth.lQm2b --grid 0:114.592:180'

Running critical command 'mktemp ABB.pot.extrapol.XXXXX'

Running subscript 'potential_extrapolate.sh --type angle ABB.pot.smooth.lQm2b ABB.pot.extrapol.CUPup' (from tags potential extrapolate) dir /opt/votca/share/votca/scripts/inverse

Extrapolate ABB.pot.smooth.lQm2b to ABB.pot.extrapol.CUPup

Running subscript 'table_extrapolate.pl --function linear --avgpoints 3 --region leftright ABB.pot.smooth.lQm2b ABB.pot.extrapol.CUPup' (from tags table extrapolate) dir /opt/votca/share/votca/scripts/inverse

Use of uninitialized value in subtraction (-) at /opt/votca/share/votca/scripts/inverse/table_extrapolate.pl line 227.

Use of uninitialized value in subtraction (-) at /opt/votca/share/votca/scripts/inverse/table_extrapolate.pl line 227.

Illegal division by zero at /opt/votca/share/votca/scripts/inverse/table_extrapolate.pl line 227.

The question is if it was really intended to be resampled with the stepsize 114.592 (see the bold line above) or it is a typo?

[Christoph Junghans]

Ok, I think I know what went wrong. The step size given in xml file
has to be in VOTCA units (radian).
And 2.0 in radian is a bit big. You want to go something like 0.034
radian = 2 degree.

Christoph

[Vitalie Botan]

Indeed, it scaled the step as well, changing to radians resolves this problem. However, there is another one I discovered today when trying to run the following command:

csg_call --options settings.xml rdf pot A-A.dist.new A_A.pot

ERROR:                                                         

source_wrapper: Could not get any script from tags 'rdf' 'pot' 

Details can be found above 

The script itself is in the usual location. Any ideas?

--Vitalie

[Christoph Junghans]

2014-12-05 7:45 GMT-07:00 Vitalie Botan <vitali...@gmail.com>:
> Indeed, it scaled the step as well, changing to radians resolves this
> problem. However, there is another one I discovered today when trying to run
> the following command:
>
> csg_call --options settings.xml rdf pot A-A.dist.new A_A.pot
>
> ERROR:
> source_wrapper: Could not get any script from tags 'rdf' 'pot'
> Details can be found above
>
> The script itself is in the usual location. Any ideas?

I am not quite sure what you want to achieve with that call.

I guess what you want is
$ csg_call dist invert

see
$ csg_call dist invert --help
and
$ csg_call -l
for all possible calls.

Christoph

[Vitalie Botan]

Hi Christoph,

It came to my attention only recently that the --min option of the command above does change the potential barrier hight. It is not clear from the description which minimum value is meant. Could you please comment what exactly does it mean:

csg_call dist invert --min 0.001

Thanks.

Best,

--Vitalie

[Christoph Junghans]

2015-05-05 18:47 GMT-06:00 Vitalie Botan <vitali...@gmail.com>:
> Hi Christoph,
>

> It came to my attention only recently that the --min option of the command
> above does change the potential barrier hight. It is not clear from the
> description which minimum value is meant. Could you please comment what
> exactly does it mean:
>
> csg_call dist invert --min 0.001

Basically, the Boltzmann inversion is U(x)=-k_B*T log(P(x)/norm(x))
with x being r (or the angle, bond length), which is undefined for
P(x)=0, so in VOTCA we only invert for P(x)>min. Because of that the
absolute barrier height is not a very good measure and relative
difference in barriers are preferred.

Christoph