derivation of potentials for dihedral
Ali Khodayari
Dears,
I am trying to derive potentials for a CG cellobiose-water system. I am using IBI. The non-bonded interactions are derived using the distributions taken from csg_stat, and I am trying to use csg_boltzmann to get the potentials for bonded interactions.
As the simulations crash at the very first step of the iteration, I tried to add the potentials one by one (first bonds, then angles, and at the end dihedrals).
An error which I face regularly is:
A tabulated bond interaction table number 1 is out of the table range: r
0.501462, between table indices 501 and 502, table length 501
I believe this is due to high forces between atoms at some point, shooting the atoms to distances higher than the table’s length.
After solving the issues with the bonds and angles, I am now trying to get a nice potential for the dihedrals in my system. Yet, the potentials for the dihedrals always look very weird (see table_d1.xvg). Without including this table, the simulations in gromacs at least run for 1000 steps (1ps), but including this crashes the simulations at the very first step. I am using the following commands for csg_boltzmann:
csg_boltzmann --cg "cellobiose.xml;water.xml" --top md_at.tpr --trj md_at.trr < boltzmann_cmds2
having Boltzmann_cmds as:
tab set T 300
tab set smooth_pot 1
tab set n 101
tab dihedral_yrry.pot *:dihedral-YRRY:*
q
The output is already looking very strange (see dihedral_yyry.pot)
Then I use the following commands to get the table:
BOND_FILE=$"dihedral_yrry.pot"
SETTINGS="convert_dihedral.xml"
# cut bad sampled regions at the boundaries
sed -e '1,33d' -e 's/$/ i/' $BOND_FILE | tac | sed -e '1,0d' | tac > "$BOND_FILE.cut"
csg_call table smooth "$BOND_FILE.cut" "$BOND_FILE.smooth"
csg_resample --in "$BOND_FILE.smooth" --out "$BOND_FILE.refined" --grid 0:0.001:3.135
csg_call table extrapolate --function quadratic --region left "$BOND_FILE.refined" "$BOND_FILE.refined1"
csg_call table extrapolate --function quadratic --region right "$BOND_FILE.refined1" "$BOND_FILE.cur"
csg_call --ia-type dihedral --ia-name dihedral --options $SETTINGS convert_potential gromacs "$BOND_FILE.cur" table_d${2}.xvg
and the .xml file looks like this:
<cg>
<bonded>
<bondtype>dihedral</bondtype>
<name>dihedral</name>
<min>0</min>
<max>3.1415</max>
<step>0.001</step>
<inverse>
<gromacs>
<table>table_d1.xvg</table>
</gromacs>
</inverse>
</bonded>
<inverse>
<gromacs>
<pot_max>1e6</pot_max>
<table_end>3.1415</table_end>
<min>0.0</min>
<max>3.1415</max>
<table_bins>0.004</table_bins>
</gromacs>
</inverse>
</cg>
I have tried to cut different edges of the potential, but they all seem not logical to me, and moreover, cause the simulation to crash at the very beginning.
Would you please kindly let me know if I am doing something wrong and/or how I can resolve this issue?
My best,
Ali
Christoph Junghans
dihedral potential conversion see:
https://github.com/votca/csg/pull/500
but in general dihedral IBI is pretty still pretty experimental in IBI
and we should most likely add a warning about that to the code.
I, personally, have not worked with dihedral potentials enough to give
any advice on this matter.
Christoph
>
>
>
> My best,
>
> Ali
>
> --
> 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/003b01d5e813%2416be71a0%24443b54e0%24%40gmail.com.
--
Christoph Junghans
Web: http://www.compphys.de
Ali Khodayari
Thank you for your response. Shall I replace the new potential_to_gromacs.sh with the one I already have?
I have another question. Could it be because of the type of dihedral I chose that the potential does not exhibit high peaks on the two ends? The thing that does not show it good to me is that I have very low energies on the two ends of the parabolic shape.
How can I adjust the amplitude of the two sides? For some of the potentials I see very large peaks of 50000 and for some only 100 after extrapolating the two sides. I also tried diffent extrapolators for my potentials.
My regards,
Ali
Christoph Junghans
>
> Dear Christoph,
>
> Thank you for your response. Shall I replace the new potential_to_gromacs.sh with the one I already have?
>
> I have another question. Could it be because of the type of dihedral I chose that the potential does not exhibit high peaks on the two ends? The thing that does not show it good to me is that I have very low energies on the two ends of the parabolic shape.
>
> How can I adjust the amplitude of the two sides? For some of the potentials I see very large peaks of 50000 and for some only 100 after extrapolating the two sides. I also tried diffent extrapolators for my potentials.
points on each side and that sometimes over/underestimates the slope,
but usually the extrapolated part of the potential doesn't really
matter much (for bonds and angles at least).
Christoph