Tabulated vs harmonic potential for a bond

[Alexander Alexander]

Dear all,
I did a
My IBI-bonded-nonbonded (after getting the nonbonded potentials via IBI-nonbonded only) crashes dues to the "...bonded interactions could not be calculated because some atoms
involved moved further apart ...". By replacing the ill-bond by a harmonic potential in topol.top file "2 3   N  0.3453  K" where N = 1,2 and K = 5000, the simulation are going well. Now what? what is the relation between this harmonic  and the table_b2.xvg which was replaced by the harmonic one? What would happen for the table_b2.xvg? Figures in attachment are the plot of table_b2.xvg.
I got the bonded tables including the table_b2.xvg by "csg_boltzmann" --> "csg_call --sloppy-tables" ---> "csg_call --ia-type bond --ia-name bond ... gromacs table_b2.xvg". And I am using the bond(angle)-*.pot.in as initial guesses for my IBI-all.

The table_b1.xvg in gromacs manual has been defined such that the the r(nm), f(r), f'(r) are the first, second and third column of table_b2.xvg where f(r) is a cubic spline function in V_b(r_ij) = k*f(r_ij), If the same definition is applied for the table_b2.xvg obtained by the above method? What is the unit of f(r)?
What is the unit of each column in table_b2.xvg in VOTCA?

And would you please explain table_b2-1-3.png figure?
Thank you very much.
Regards,
Alex

[Christoph Junghans]

On Wed, May 2, 2018 at 1:38 PM, Alexander Alexander
<alexand...@gmail.com> wrote:
> Dear all,
> I did a
> My IBI-bonded-nonbonded (after getting the nonbonded potentials via
> IBI-nonbonded only) crashes dues to the "...bonded interactions could not be
> calculated because some atoms
> involved moved further apart ...". By replacing the ill-bond by a harmonic
> potential in topol.top file "2 3 N 0.3453 K" where N = 1,2 and K = 5000,
> the simulation are going well. Now what? what is the relation between this
> harmonic and the table_b2.xvg which was replaced by the harmonic one? What
> would happen for the table_b2.xvg?

Sorry, I don't understand your question! Can you rephrase?

> Figures in attachment are the plot of
> table_b2.xvg.
> I got the bonded tables including the table_b2.xvg by "csg_boltzmann" -->
> "csg_call --sloppy-tables" ---> "csg_call --ia-type bond --ia-name bond ...
> gromacs table_b2.xvg". And I am using the bond(angle)-*.pot.in as initial
> guesses for my IBI-all.
>
> The table_b1.xvg in gromacs manual has been defined such that the the r(nm),
> f(r), f'(r) are the first, second and third column of table_b2.xvg where
> f(r) is a cubic spline function in V_b(r_ij) = k*f(r_ij), If the same
> definition is applied for the table_b2.xvg obtained by the above method?

To be clear, the gromacs table format is x, f(x), -f'(x) (see gromacs
manual 4.2.14)

I am not sure what your question is! Gromacs uses cubic spline
internal to interpolate the table.
However, if you set cg.inverse.gromacs.table_bins in your xml file
small enough (you had something like 0.002), then there isn't much to
interpolate.

> What is the unit of f(r)?

KJ/mol (gromacs energy units)

> What is the unit of each column in table_b2.xvg in VOTCA?

KJ/mol/nm (gromacs force units)

>
> And would you please explain table_b2-1-3.png figure?

Not sure what there is to explain..the code on how the third column is
calculated is here:
https://github.com/votca/csg/blob/master/share/scripts/inverse/table_to_xvg.pl#L95

Christoph

> Thank you very much.
> Regards,
> Alex
>

> --
> 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 https://groups.google.com/group/votca.
> For more options, visit https://groups.google.com/d/optout.



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

[Alexander Alexander]

Hi Christoph,
Thanks.

On Wednesday, May 2, 2018 at 5:13:13 PM UTC-4, Christoph Junghans wrote:

On Wed, May 2, 2018 at 1:38 PM, Alexander Alexander
<alexand...@gmail.com> wrote:
> Dear all,
> I did a
> My IBI-bonded-nonbonded (after getting the nonbonded potentials via
> IBI-nonbonded only) crashes dues to the "...bonded interactions could not be
> calculated because some atoms
> involved moved further apart ...". By replacing the ill-bond by a harmonic
> potential in topol.top file "2 3   N  0.3453  K" where N = 1,2 and K = 5000,
> the simulation are going well. Now what? what is the relation between this
> harmonic  and the table_b2.xvg which was replaced by the harmonic one? What
> would happen for the table_b2.xvg?
Sorry, I don't understand your question! Can you rephrase?

Let's say the calculation crashes when table_b2.xvg is used but it works fine with "2 3   2  0.3453  1000" instead, then my question is that "should I continue the preregistration using "2 3   2  0.3453  1000" and forget about the table_b2.xvg? If so, what would be the final right table_b2.xvg at the end? If not so, how I can take benefit from the knowledge of  "2 3   2  0.3453  1000" to improve the table_b2.xvg so that it would work as well?

 

> Figures in attachment are the plot of
> table_b2.xvg.
> I got the bonded tables including the table_b2.xvg by "csg_boltzmann" -->
> "csg_call --sloppy-tables" ---> "csg_call --ia-type bond --ia-name bond ...
> gromacs table_b2.xvg". And I am using the bond(angle)-*.pot.in as initial
> guesses for my IBI-all.
>
> The table_b1.xvg in gromacs manual has been defined such that the the r(nm),
> f(r), f'(r) are the first, second and third column of table_b2.xvg where
> f(r) is a cubic spline function in V_b(r_ij) = k*f(r_ij), If the same
> definition is applied for the table_b2.xvg obtained by the above method?
To be clear, the gromacs table format is x, f(x), -f'(x) (see gromacs
manual 4.2.14)

I am not sure what your question is! Gromacs uses cubic spline
internal to interpolate  the table.
However, if you set cg.inverse.gromacs.table_bins in your xml file
small enough (you had something like 0.002), then there isn't much to
interpolate.


> What is the unit of f(r)?
KJ/mol (gromacs energy units)

> What is the unit of each column in table_b2.xvg in VOTCA?
KJ/mol/nm (gromacs force units)

Just to sum up the table format, can you please confirm me that the table_b.xvg in VOTCA has the following format:
first column : Distance with unit nm
Second column : Energy with unit KJ/mol
Thirst column: Force with unit KJ/mol/nm 

>
> And would you please explain table_b2-1-3.png figure?
Not sure what there is to explain..the code on how the third column is
calculated is here:
https://github.com/votca/csg/blob/master/share/scripts/inverse/table_to_xvg.pl#L95

Best regards,
Alex

[Christoph Junghans]

On Wed, May 2, 2018 at 3:37 PM, Alexander Alexander

You know generate a table with k*(r-r_0) (with r_0=0.3453 andk=1000,
clip out the middle part and replace it with the structured part from
BI.

For bonds (not angles) in gromacs that is the format!

Also for a different simulation backend, e.g. lammps this will be
different of course, too!

Christoph

[Alexander Alexander]

Thanks.

OK, I see.
How about this instead:
The reason that the calculation crashes with table_b2.xvg is that the table_b2.xvg represents a very low K (spring constant) so that the two atoms are getting apart. I was reading somewhere in the mailing list that the last number in "2 3   8  2  1  ; 1:bond-H1-L:1" is a kind of factor for the third column of the table_b2.xvg, so, don't you think it would be possible to increase this factor to 2 or 3 (which means multiplying the third column) until the calculation goes fine?

For angle the first column is degree probably.

Also for a different simulation backend, e.g. lammps this will be
different of course, too!

I have one more question which your comment on would be highly appreciated:
In another IBI-bonded-nonbonded parameterization the "mpirun -np 64 gmx_mpi mdrun" work fine till step_021 of parameterization where it crashes because it can not find any domain decomposition for 64 ranks that is compatible with the system size. Even it works in the minimization simulation in step_021.
I have already read most of the discussions about the error of domain decomposition, and I have tried different options like using different ranks (64, 32, 16, 8, 40, 36) and also -rdd ... or -dds. However, no success yet. I am totally confused with these parameters and I can not solve the issue. So, below I have shared the md.log files for the step_20 (which works fine) and step_21 (minimization which works fine) and step_21 (which crashes), and I would be so appreciated if one could help me overcome the issue.
About the machine: Running on 4 nodes with total 64 cores, 64 logical cores, Cores per node: 16.

md.log-20
https://drive.google.com/open?id=1ZwYggDI1Hy_eYY8_nT8I50t91Wla0Obr

md.log-21-mini
https://drive.google.com/open?id=1oYE4A6fC1L8u_vx4sR8RIq3EWiG9-7YB

md.log-21
https://drive.google.com/open?id=1ff_QJrdGvLFgwGU9xMY-IcqfxzF2jkAO

[Christoph Junghans]

On Wed, May 2, 2018 at 4:37 PM, Alexander Alexander

Yeah, that sounds reasonable! Increasing the K helps usually helps
with the atoms getting too far apart.

If you want to scale column 3 you will also have to modify column 2 accordingly.
Gromacs won't allow mismatching columns.

Sorry, that really is a question for the gromacs user mailing list.
(gromacs.or...@maillist.sys.kth.se)

[Alexander Alexander]

Hi Christoph,

I tried to scale both of the column 2 and 3 of table_b2.xvg using a lot of scale-factor, but they did not work and the missing bond ... error is still here., I guess finding the right factor is not easy.
So, I tried your suggestion above, I mean k*(r-r_0) (with r_0=0.3453 andk=1000. T do so, I first created a bond.pot.ib potential including 3 columns as following:
#-----------------
bond.pot.in
#column1          #column2                          # column3
r                (1000/2)*(r-0.3453)^2         -(1000)*(r-0.3453)
#-------------------
And added the letter "i" to the end of each row.
Then I tried to converted the above bond.pot.in to table_b2.xvg using "csg_call --ia-type bond --ia-name bond-H1-L --options bond-H1-L.xml convert_potential gromacs ...". But then all of the second and third columns of the new table_b2.xvg was just nan -nan.
Would you please let me know if the procedure is right and where I am doing wrong?

Thank you.
Alex

[Christoph Junghans]

On Thu, May 3, 2018 at 12:44 PM, Alexander Alexander

Not sure, can you post your input files? (pot.in and xml)

Christoph

[Alex]

Thanks.

Please find attached the files pot.in =  Harmonic-bond-H1-L.pot.in

The two other files are the potential associated with this bond obtained by normal BI crashing the calculation.

 I simply calculate the harmonic files using these two awk command:

awk -v CONVFMT=%.7g '{$3 = (($1 - 0.34530000)**2)*(1000/2)} 1'

awk -v CONVFMT=%.7g '{$4 = ($1 - 0.34530000)*(-1000)} 1'

Regards,

Alex

> email to votca+unsubscribe@googlegroups.com.

> To post to this group, send email to vo...@googlegroups.com.
> Visit this group at https://groups.google.com/group/votca.
> For more options, visit https://groups.google.com/d/optout.



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

--

You received this message because you are subscribed to a topic in the Google Groups "votca" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/votca/7Hh-gxZEpmI/unsubscribe.
To unsubscribe from this group and all its topics, send an email to votca+unsubscribe@googlegroups.com.

[Christoph Junghans]

Ok, I think there is some confusion about table format.

The Gromacs (xvg) table format is:

first column : Distance with unit nm
Second column : Energy with unit KJ/mol
Thirst column: Force with unit KJ/mol/nm

However the VOTCA table format is just
first column : Distance in simulation backend units
Second column : Energy in simulation backend units
Third column: flag
(see section 3.6 of the VOTCA manual)

>> > Then I tried to converted the above bond.pot.in to table_b2.xvg using
>> > "csg_call --ia-type bond --ia-name bond-H1-L --options bond-H1-L.xml
>> > convert_potential gromacs ...". But then all of the second and third
>> > columns
>> > of the new table_b2.xvg was just nan -nan.
>> > Would you please let me know if the procedure is right and where I am
>> > doing
>> > wrong?
>> Not sure, can you post your input files? (pot.in and xml)
>
> Thanks.
> Please find attached the files pot.in = Harmonic-bond-H1-L.pot.in
> The two other files are the potential associated with this bond obtained by
> normal BI crashing the calculation.
> I simply calculate the harmonic files using these two awk command:
>
> awk -v CONVFMT=%.7g '{$3 = (($1 - 0.34530000)**2)*(1000/2)} 1'
> awk -v CONVFMT=%.7g '{$4 = ($1 - 0.34530000)*(-1000)} 1'

Looking at Harmonic-bond-H1-L.pot.ib, something is screwed up.
You have duplicated r values.
$ head Harmonic-bond-H1-L.pot.ib
0.3432 0.002205 2.1
0.3432 0.002205 2.1

Something like:
$ awk 'BEGIN{for(r=0.3425;r<0.3481;r+=0.0005){print
r,((r-0.34530000)**2)*(1000/2)},"i";}' > Harmonic-bond-H1-L.pot.ib
$ csg_call --ia-type bond --ia-name bond-H1-L --options bond-H1-L.xml
convert_potential gromacs Harmonic-bond-H1-L.pot.ib table_b2.xvg
should work!

Christoph

>> You received this message because you are subscribed to a topic in the
>> Google Groups "votca" group.
>> To unsubscribe from this topic, visit
>> https://groups.google.com/d/topic/votca/7Hh-gxZEpmI/unsubscribe.
>> To unsubscribe from this group and all its topics, send an email to

[Alexander Alexander]

Probably you mena the *.pot.ib or *.pot.in we get from csg_boltzmann. Otherwise after converting them to the gromacs format by by csg_call... .. convert_potential gromacs... they have the gromacs units and formats. 

Wow, it works now, thank you.
Just one issue; in "2 3   2  0.3453  K" of harmonic potential , different K from 500 to 50000 could work in this system, and then each K would give different pot.ib and table_b2.xvg of course with (r-0.34530000)**2)*(1000/2), so, which one (K) is good to choose?

One more question:
After a IBI-nonbonded for a system, I am doing now IBI-bonded-nonbonded for that system which has two bond types and one angle. They started with a promising fit with their dist.tgt's but afterward they do not have any sensible improvement after 115 step now, even the angle.dist.new is loosing gradually its agreement with angle.dist.tgt. I was wondering what might be the reason? Below is part of my settings.xml file: The <scale>1.5</scale> in the first bond I choose just because I though it might help improvement, but no considerable improvement.

     <name>bond-D-E</name>
     <min>0.200</min>
     <max>0.500</max>
     <step>0.0005</step>
           <inverse>
                <target>bond-D-E.dist.tgt</target>
                <do_potential>1 0 0</do_potential>
                <post_update>smooth</post_update>
                <post_update_options>
                <scale>1.5</scale>
                </post_update_options>
                <post_add>convergence</post_add>
                <gromacs>
                <table>table_b1.xvg</table>
                </gromacs>
           </inverse>
  </bonded>
<!--  %%% -->
  <bonded>
     <name>angle-D-E-F</name>
     <min>2.0</min>
     <max>3.0</max>
     <step>0.001</step>
          <inverse>
                <target>angle-D-E-F.dist.tgt</target>
                <do_potential>0 0 1</do_potential>
                <post_update>smooth</post_update>
                <post_update_options>
                <scale>0.1</scale>
                </post_update_options>
                <post_add>convergence</post_add>
                <gromacs>
                <table>table_a1.xvg</table>
                </gromacs>
         </inverse>
  </bonded>

Thank you.
Regards,
Alex

[Alexander Alexander]

Hi Christoph,

Any idea, please?

Additionally, why using the "2 3   2  0.3453  K" in the harmonic format works fine in a calculation, but it does not work soemtimes when I convert the "0.3453  K" to the table format following the remedy above?

Best regards,
Alex

[Christoph Junghans]

On Sat, May 5, 2018 at 5:46 AM, Alexander Alexander
<alexand...@gmail.com> wrote:
> Hi Christoph,
>
> Any idea, please?
I am out of ideas for now!

Have another look here:
http://www.gromacs.org/Documentation/How-tos/Tabulated_Potentials#Constructing_the_Table(s)
If you are doing the table right, otherwise this sounds like a bug in
gromacs to me!

Christoph

[Jakub Krajniak]

Hi Alexander, 

In the very first e-mail, you have shown us the plot of table_.xvg, I suppose this is the forces column. Now you can see that the force is 0 before and after the region of interest. 

In fact, it does not matter how big K you put in topol.top for bond type 8, K*0 = 0. I think that's why the atoms fly away. What I usually did is to extrapolate the tables outside the interested regions.

So please, check your tables, the forces column. If it contains a sequence of 0 for low values of r then you have a problem here.

I use following script to do the conversion from a pure potential file (distance, energy, flag) to GROMACS bond tables:

https://github.com/md-lab-tools/lab-tools/blob/master/scripts/votca/bond_convert.sh

to convert .pot files to GROMACS format with the setting file:

https://github.com/md-lab-tools/lab-tools/blob/master/scripts/votca/convert_bond.xml

Also, I have noticed that you use interchangeably bond types 1 and 2. Keep in mind that only bond type 1 is harmonic and bond type 2 is quadratic (so-called GROMOS bonds).

Best,

Jakub

[Christoph Junghans]

On Sat, May 5, 2018 at 09:03 Jakub Krajniak <jkra...@gmail.com> wrote:

Hi Alexander, 

In the very first e-mail, you have shown us the plot of table_.xvg, I suppose this is the forces column. Now you can see that the force is 0 before and after the region of interest. 

In fact, it does not matter how big K you put in topol.top for bond type 8, K*0 = 0. I think that's why the atoms fly away. What I usually did is to extrapolate the tables outside the interested regions.

So please, check your tables, the forces column. If it contains a sequence of 0 for low values of r then you have a problem here.

I use following script to do the conversion from a pure potential file (distance, energy, flag) to GROMACS bond tables:

https://github.com/md-lab-tools/lab-tools/blob/master/scripts/votca/bond_convert.sh

to convert .pot files to GROMACS format with the setting file:

https://github.com/md-lab-tools/lab-tools/blob/master/scripts/votca/convert_bond.xml

Also, I have noticed that you use interchangeably bond types 1 and 2. Keep in mind that only bond type 1 is harmonic and bond type 2 is quadratic (so-called GROMOS bonds).

Best,

Jakub

Thanks Jakub!

Christoph 

[Alexander Alexander]

Hi,

On Saturday, May 5, 2018 at 11:03:52 AM UTC-4, Jakub Krajniak wrote:

Hi Alexander, 

In the very first e-mail, you have shown us the plot of table_.xvg, I suppose this is the forces column. Now you can see that the force is 0 before and after the region of interest. 

In fact, it does not matter how big K you put in topol.top for bond type 8, K*0 = 0. I think that's why the atoms fly away. What I usually did is to extrapolate the tables outside the interested regions.

So please, check your tables, the forces column. If it contains a sequence of 0 for low values of r then you have a problem here.

I use following script to do the conversion from a pure potential file (distance, energy, flag) to GROMACS bond tables:

https://github.com/md-lab-tools/lab-tools/blob/master/scripts/votca/bond_convert.sh

to convert .pot files to GROMACS format with the setting file:

https://github.com/md-lab-tools/lab-tools/blob/master/scripts/votca/convert_bond.xml

Thanks Jakub for the clarification.
Actually my my potential or table do not have a sequence of 0, it is 0 only in the first point.
Anyway, I started the parametrization from scratch and using your script I extended all of the potentials and used those new potentials as *.pot.in in my mean directory. However, the calculation was crashed in stap_006, even earlier that my old calculation. The attached are the potentials of the bond which cause the problem. And below is the setting of it:
  <bonded>
     <name>bond-H1-L</name>
     <min>0.3425</min>
     <max>0.348</max>
     <step>0.0005</step>
           <inverse>
                <target>bond-H1-L.dist.tgt</target>
                <do_potential> 0 1 1</do_potential>

                <post_update>smooth</post_update>
                <post_update_options>
                <scale>0.1</scale>
                </post_update_options>
                <post_add>convergence</post_add>
                <gromacs>

                <table>table_b2.xvg</table>
                </gromacs>
           </inverse>
  </bonded>
The bond-H1-L.pot.in is 0 to 0.6 nm but due to the setting they will be tuned to the limit defined by the setting, and that might be the problem, would you please let me know what do you think?
Regards,
Alex

Although the settings.xml file limit the potential in the small region of interest.