Energy divergence when performing coarse-grained simulation

[Satyen Dhamankar]

Hello,

I am performing a coarse-grained simulation of benzene, where I am converting C6H6 into two beads. Where each bead is essentially C3H3, and they are connected by a bond. 

I want to perform IBI on this coarse-grained simulation using votca and gmx, but I am having some problems. 

The command I am running is:
csg_inverse --options settings.xml

I am getting an energy divergence (blows up to infinity) when votca runs mdrun. This is the error I see in my inverse.log file:

Program: gmx mdrun, version 2019.6
Source file: src/gromacs/mdlib/sim_util.cpp (line 752)
MPI rank: 1 (out of 6)
Fatal error:
Step 0: The total potential energy is inf, which is not finite. The LJ and
electrostatic contributions to the energy are inf and 0, respectively. A
non-finite potential energy can be caused by overlapping interactions in
bonded interactions or very large or Nan coordinate values. Usually this is
caused by a badly- or non-equilibrated initial configuration, incorrect
interactions or parameters in the topology.
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors

I am running a minimization script, but I don't see why energy is blowing up. 

I also am unsure of what the table_b1.xvg and table_a1.xvg files are for. I naively got these files from the hexane tutorial. What are these files for? How do I make them?

I have attached my grommp.mdp, conf.gro, settings.xml, topol.top and tables files.  

Any advice you have would be appreciated! 

[Christoph Junghans]

These two tables are for the bonded interactions, they are just
boltzmann-inverted tables from the bonded distributions (do that by
hand or use csg_boltzmann).
You could also use harmonic potentials as a start to make your life easier.

Christoph

>
> I have attached my grommp.mdp, conf.gro, settings.xml, topol.top and tables files.
>
> Any advice you have would be appreciated!
>
>

> --
> Join us on Slack: https://join.slack.com/t/votca/signup
> ---
> 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 view this discussion on the web visit https://groups.google.com/d/msgid/votca/956dbed6-f6a1-4269-866c-a4a8c5a6428cn%40googlegroups.com.



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

[Satyen Dhamankar]

Thank you for your advice, Christoph. 

How do you recommend I use harmonic potentials for a CG simulation? I have the bond.dist.new histogram after I did my mapping, I can use the modal bond length I obtained from that. However, what can I use as my spring constant? 

Using csg_boltzmann on my bond distribution, I obtained the table at 298 K and set bond, however, the table that was created started at 0.19 nm. My question is, when performing a boltzmann inversion on the bond distribution, how do we treat the bins with a height of 0? should the number at that point simply be inf? 

This is my mapping file, just to be complete:

<cg_molecule>

<name>BNZ</name> <!-- the name of the molecule after mapping--> 

<ident>UNK</ident> <!-- name of the actual atom --> 

<topology>

<cg_beads>

<cg_bead>

<name>B1</name>

<type>CG</type>

<symmetry>1</symmetry>

<mapping>A</mapping> 

<beads>

1:UNK:C00 1:UNK:C01 1:UNK:C02 1:UNK:H06 1:UNK:H07 1:UNK:H08 

</beads> 

</cg_bead>

<cg_bead> 

<name>B2</name>

<type>CG</type>

<symmetry>1</symmetry> 

<mapping>A</mapping> 

<beads>

1:UNK:C03 1:UNK:C04 1:UNK:C05 1:UNK:H09 1:UNK:H0A 1:UNK:H0B

</beads>

</cg_bead>

</cg_beads>

<cg_bonded> 

<bond>

<name>bond</name> 

<beads>

B1 B2

</beads>

</bond>

</cg_bonded>

</topology>

<maps>

<map>

<name>A</name>

<weights>12.0110 12.0110 12.0110 1.0080 1.0080 1.0080</weights>

</map>

</maps>

</cg_molecule>

I used this file as my mapping file when I executed:

csg_boltzmann --top topol.tpr --trj trajectory.trj --cg bezene.xml

tab set T 298

tab set scale bond 

tab table_b1.xvg *:bond:*

[Christoph Junghans]

On Mon, Apr 5, 2021 at 22:56 Satyen Dhamankar <saty...@gmail.com> wrote:

Thank you for your advice, Christoph. 

How do you recommend I use harmonic potentials for a CG simulation? I have the bond.dist.new histogram after I did my mapping, I can use the modal bond length I obtained from that. However, what can I use as my spring constant? 

you can do a boltzmann inversion and then do a quadratic fit.

Using csg_boltzmann on my bond distribution, I obtained the table at 298 K and set bond, however, the table that was created started at 0.19 nm. My question is, when performing a boltzmann inversion on the bond distribution, how do we treat the bins with a height of 0? should the number at that point simply be inf?

You will have to do a little bit of post-processing: https://www.votca.org/csg/preparing.html#post-processing-of-the-potential

To view this discussion on the web visit https://groups.google.com/d/msgid/votca/5b85becb-b3ce-414e-8409-d6f2399e7a81n%40googlegroups.com.

[Satyen Dhamankar]

Hello,

The simulation I am running is of a molecule comprising of two identical beads held together by a bond. 

I ran the following on my coarse-grained simulation:

csg_inverse --options settings.xml --debug

I have been working on this, and I have been receiving this message:

Appending to existing logfile inverse.log
We are doing Method: ibi
Prepare (dir step_000)
Using initial guess from dist bond.dist.tgt for bond
#############################################################################################################################################
# #
# ERROR: #
# do_external: subscript #
# /home/satyend/software/spack/opt/spack/linux-rhel7-broadwell/gcc-9.3.1/votca-csg-2021-vh2kxfbymvudsn3en3wk5xxnie4j73nq/share/votca/scripts/inverse/dist_boltzmann_invert.pl --type bond --kbT 2.47772 --min 1e-10 bond.dist.tgt bond.pot.new.raw.cQ7 #
# (from tags dist invert) failed #
# For details see the logfile /tigress/satyend/benzene_simulations/opls-aa/r_cutoff-1.5-org/VOTCA/TwoBead_CG/ibi/inverse.log #
##############################################################################################################################################
Terminated

What could the cause of such an error message be? 

[Satyen Dhamankar]

This is also something that I saw... 
123: die "All data points from file '$infile' are invalid after Boltzmann inversion, please check if your distribution is a valid dist.\n" if ($valid_i==-1);
126: my $first=undef;
127: for (my $i=$valid_i;$i>=0;$i--){
128: if ($flag[$i] eq "u") {
134: $first=0 unless defined $first;
137: my $last=undef;
138: for (my $i=$valid_i;$i<=$#pot;$i++){
139: if ($flag[$i] eq "u") {
139: if ($flag[$i] eq "u") {
139: if ($flag[$i] eq "u") {
139: if ($flag[$i] eq "u") {
139: if ($flag[$i] eq "u") {
139: if ($flag[$i] eq "u") {
139: if ($flag[$i] eq "u") {
145: $last=$#pot unless defined $last;
146: my $n=10;
147: my $valid=$last-$first+1;
148: die "Only $valid points are valid after Boltzmann inversion from file '$infile', please check if your distribution is a valid dist.\n" if ($valid < $n);

Does this imply that I have to run post-processing? My question is, I see the same zeros in the hexane bond.dist.tgt file, but I did not receive an error when I ran that simulation: why is that?

[Christoph Junghans]

On Wed, Apr 7, 2021 at 20:22 Satyen Dhamankar <saty...@gmail.com> wrote:

This is also something that I saw... 
123: die "All data points from file '$infile' are invalid after Boltzmann inversion, please check if your distribution is a valid dist.\n" if ($valid_i==-1);
126: my $first=undef;
127: for (my $i=$valid_i;$i>=0;$i--){
128: if ($flag[$i] eq "u") {
134: $first=0 unless defined $first;
137: my $last=undef;
138: for (my $i=$valid_i;$i<=$#pot;$i++){
139: if ($flag[$i] eq "u") {
139: if ($flag[$i] eq "u") {
139: if ($flag[$i] eq "u") {
139: if ($flag[$i] eq "u") {
139: if ($flag[$i] eq "u") {
139: if ($flag[$i] eq "u") {
139: if ($flag[$i] eq "u") {
145: $last=$#pot unless defined $last;
146: my $n=10;
147: my $valid=$last-$first+1;
148: die "Only $valid points are valid after Boltzmann inversion from file '$infile', please check if your distribution is a valid dist.\n" if ($valid < $n);

Does this imply that I have to run post-processing?

correct!

My question is, I see the same zeros in the hexane bond.dist.tgt file, but I did not receive an error when I ran that simulation: why is that?

you also need to check the flags (3rd column) in the Distribution file as well , „u“ means undefined, you want an „i“ there for in range.

Christoph 

To view this discussion on the web visit https://groups.google.com/d/msgid/votca/51ef0700-6df2-43d9-8f62-2a5d6f196e76n%40googlegroups.com.

[Satyen Dhamankar]

I apologize for the incessant questioning... 

But this is the procedure i have followed: 

I ran a all-atom simulation on gromacs, obtained all the necessary file (topol.trr, conf.gro, ener.edr, traj.trr). 

I wrote a mapping file to turn benzene into a two-bead structure. 
Once I got these files, I ran the following command:
csg_boltzmann --top 600-c6h6-MolDynamics.tpr --trj 600-c6h6-MolDynamics.trr --cg benzene.xml < boltzmann_cmds

where boltzmann_cmds is:

hist bond.dist.ib *bond*
q

I got bond.dist.ibi, bond.dist.new from this. I also received A-A.dist.new from:

csg_stat --top 600-c6h6-MolDynamics.tpr --trj 600-c6h6-MolDynamics.trr --cg benzene.xml --options settings.xml --nt 6 --begin 200

My bond.dist.new has no u's in the third column. However, it has values only from 0.192 to 0.204 - so a very small range. 

This is what I am using as my bond.dist.tgt in my IBI. Everything has the "i" flag in the third column. 

I am getting 3 perl_debug.xxx files in my step_000 folder. The last one had the information I have given in the above message. 

As a matter of fact, I ran the same procedure with the hexane file, and it checked out just fine... 

I am now wondering if there is something amiss with my settings.xml file, but the errors seem to be pointing squarely at bond.dist.tgt file... 

Could you explain what I have to do preprocess the data? 

I did what the page you mentioned asked me to do: 

csg_resample --in bond.dist.new --out bond.dist.newer --grid 0:0.001:1

I have written out the debug files here: https://justpaste.it/3q0xm

Any advice you have would be appreciated!! 

[Christoph Junghans]

On Wed, Apr 7, 2021 at 11:23 PM Satyen Dhamankar <saty...@gmail.com> wrote:
>
> I apologize for the incessant questioning...
> But this is the procedure i have followed:
> I ran a all-atom simulation on gromacs, obtained all the necessary file (topol.trr, conf.gro, ener.edr, traj.trr).
> I wrote a mapping file to turn benzene into a two-bead structure.
> Once I got these files, I ran the following command:
> csg_boltzmann --top 600-c6h6-MolDynamics.tpr --trj 600-c6h6-MolDynamics.trr --cg benzene.xml < boltzmann_cmds
> where boltzmann_cmds is:
>
> hist bond.dist.ib *bond*
> q
>
> I got bond.dist.ibi, bond.dist.new from this. I also received A-A.dist.new from:
> csg_stat --top 600-c6h6-MolDynamics.tpr --trj 600-c6h6-MolDynamics.trr --cg benzene.xml --options settings.xml --nt 6 --begin 200
>
> My bond.dist.new has no u's in the third column. However, it has values only from 0.192 to 0.204 - so a very small range.
> This is what I am using as my bond.dist.tgt in my IBI. Everything has the "i" flag in the third column.
>
> I am getting 3 perl_debug.xxx files in my step_000 folder. The last one had the information I have given in the above message.
>
> As a matter of fact, I ran the same procedure with the hexane file, and it checked out just fine...
>
> I am now wondering if there is something amiss with my settings.xml file, but the errors seem to be pointing squarely at bond.dist.tgt file...
> Could you explain what I have to do preprocess the data?
> I did what the page you mentioned asked me to do:
>
> csg_resample --in bond.dist.new --out bond.dist.newer --grid 0:0.001:1
>
> I have written out the debug files here: https://justpaste.it/3q0xm
>
> Any advice you have would be appreciated!!

Ok the perl files aren't super useful for find issue in the input
files, they are more for then you need to add a script to VOTCA.
Here is what I would try;
- also make
- look at inverse.log, above the big error banner there is usually an
error message from the script
- try to run the script that fails by hand, in your case, tags were
"dist invert", so try
$ csg_call --options settings.xml --ia-name bond --ia-type bond dist
invert --type bond --kbT 2.47772 --min 1e-10 bond.dist.tgt
bond.pot.new.raw.cQ7
- look at the input file (you already did that), e.g. bond.dist.tgt to
see if the range is right.

> To view this discussion on the web visit https://groups.google.com/d/msgid/votca/7dc60b54-2bb5-40a9-ac92-cd7db74696bbn%40googlegroups.com.

[Christoph Junghans]

On Thu, Apr 8, 2021 at 07:52 Christoph Junghans <jung...@votca.org> wrote:

On Wed, Apr 7, 2021 at 11:23 PM Satyen Dhamankar <saty...@gmail.com> wrote:
>
> I apologize for the incessant questioning...
> But this is the procedure i have followed:
> I ran a all-atom simulation on gromacs, obtained all the necessary file (topol.trr, conf.gro, ener.edr, traj.trr).
> I wrote a mapping file to turn benzene into a two-bead structure.
> Once I got these files, I ran the following command:
> csg_boltzmann --top 600-c6h6-MolDynamics.tpr --trj 600-c6h6-MolDynamics.trr --cg benzene.xml < boltzmann_cmds
> where boltzmann_cmds is:
>
> hist bond.dist.ib *bond*
> q
>
> I got bond.dist.ibi, bond.dist.new from this. I also received A-A.dist.new from:
> csg_stat --top 600-c6h6-MolDynamics.tpr --trj 600-c6h6-MolDynamics.trr --cg benzene.xml --options settings.xml --nt 6 --begin 200
>
> My bond.dist.new has no u's in the third column. However, it has values only from 0.192 to 0.204 - so a very small range.

Also make sure the range in your settings.xml is the same.

[Satyen Dhamankar]

Thank you Christoph, I have figured out the problem! 

Sincerely,

Satyen Dhamankar

To view this discussion on the web visit https://groups.google.com/d/msgid/votca/CAHG27e5dHtWcyVi%3DLm59UaRHkAhTvnC1yA17TQcRicmdj08HXA%40mail.gmail.com.

[Christoph Junghans]

On Thu, Apr 8, 2021 at 9:43 AM Satyen Dhamankar <saty...@gmail.com> wrote:
>
> Thank you Christoph, I have figured out the problem!

Nice. Could you please post your solution here for others to learn?

Thanks,

Christoph

> To view this discussion on the web visit https://groups.google.com/d/msgid/votca/CAFMHqF9wbuaz-SJWNgnG8NgmPdput4LuKEv_7QttHU_kEqVO-g%40mail.gmail.com.

[Satyen Dhamankar]

The problem was kind of silly. 

My settings.xml file had the inappropriate step size for my target bonds. 

so I changed my settings file to 

<bonded>
  <name>bond</name>
  <min>0.192</min>
  <max>0.204</max>
  <step>0.001</step>

.

.

.

</bonded>
My step was 0.002, instead of 0.001. Once I made that change, the program ran, and I got acceptable results. 

[Christoph Junghans]

On Thu, Apr 8, 2021 at 7:02 PM Satyen Dhamankar <saty...@gmail.com> wrote:
>
> The problem was kind of silly.
> My settings.xml file had the inappropriate step size for my target bonds.
> so I changed my settings file to
> <bonded>
> <name>bond</name>
> <min>0.192</min>
> <max>0.204</max>
> <step>0.001</step>
> .
> .
> .
> </bonded>

Ah ok, thanks for letting us know.

Christoph

> To view this discussion on the web visit https://groups.google.com/d/msgid/votca/82e3b11c-756b-43c1-820b-20e14b2257fan%40googlegroups.com.