table_*.xvg in the propane tutorial
Tinashe
tutorial.
I obtained the bond potentials using csg_boltzmann:
The following are parts of the "bond.pot.ib" and "bond.xml" files:
"bond.pot.ib"
0.152097 21.2962 5895.79
0.152405 19.4769 2445.97
...
0.165675 0.0496598 149.931
0.165983 0.0142966 80.4637
0.166292 0 17.5345
0.1666 0.00347479 -44.5805
0.166909 0.0275138 -110.011
...
0.180178 20.7073 -2300.42
0.180487 21.1912 -3816
"bond.xml"
<cg>
<inverse>
<program>gromacs</program>
<gromacs>
<pot_max>1e8</pot_max>
<table_end>3</table_end>
<table_bins>0.002</table_bins>
</gromacs>
</inverse>
</cg>
I have done the following command to obtain the table_*.xvg file:
csg_call --options bond.xml --ia-type bonded convert_potential gromacs
bond.pot.ib table_b1.xvg
and read previous suggestions to Sergio for example concerning ibi
(https://groups.google.com/forum/#!msg/votca/OSaBKQTR7C0/wLAXmbGxw74J)
but I still can't get my table_b1.xvg to match the one in the
tutorial. I would really appreciate any help to resolve this issue.
Whilst the minimum point coincides, the right and left sides of the
potentials are very steep.
How is it that the tutorial file intercepts the vertical axis at
6677.89???
0.000000 6677.890000 0
0.001000 6595.820000 0
0.002000 6514.250000 0
I have tried to modify the bin size in the bond.xml and changed the
pot_max value but my table_b1.xvg file is still different??? I surely
have to be doing something wrong, please help me.
Thanks.
Victor Ruehle
it might be the way how the unsampled areas are extrapolated. I think
the standard is linear interpolation (but i might be wrong on that).
You can always change that by using
1) csg_resample --in xxx --out yyy --grid 0:0.001:3 (choose proper
grid for your case)
2) csg_call table extrapolate --function quadtratic infile outfile
function can be linear,exponetnial, quadratic, ... check the help
3) convert potential to gromacs.
However, I don't think that will be an issue since the potential
anyways is very large and should never be sampled. In any case, plot
the potential and looks whether it looks sensible. If it has an easy
functional form, e.g. harmonic, I would recommend fitting a functional
form and use these parameters for gromacs input.
Best,
Victor
2012/2/13 Tinashe <tinash...@gmail.com>:
> --
> You received this message because you are subscribed to the Google Groups "votca" group.
> To post to this group, send email to vo...@googlegroups.com.
> To unsubscribe from this group, send email to votca+un...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/votca?hl=en.
>
Christoph Junghans
I had another look at the tutorials and updated the fmatch.xml in
propane/atomistic.
It does change much, but it make it easier to come to the table in propane/ibi.
So here is what I have done to get table_b1.xvg
$csg_stat --top topol.tpr --trj traj.trr --cg propane.xml --nt 5
--options fmatch.xml
taken from run.sh, it creates bond.dist.new, which has only points,
p>0 in the region 0.154 to 0.179.
$csg_resample --in bond.dist.new --out bond.dist.resample --grid
0.154:0.001:0.179
This crops away the unimportant region
$awk -v kbt=1.6629 '{print $1, -kbt*log($2/$1/$1)}' bond.dist.resample
> bond.pot.ib
(kbt=1.6629 as the propane example is done at T=200K) It does the
Boltzmann inversion.
$csg_call --options bond.xml --ia-type bond convert_potential gromacs
bond.pot.ib table_b1.xvg
If you plot table_b1.xvg in the region from 0.154 to 0.179 vs the
table_b1.xvg in ibi you will see that they look pretty similar.
Note1: the important region is the minimum + 2kbt.
Note2: the table_b1.xvg in ibi folder is a hand-made table done in 2008.
The rest of the deviation are due to technical things, which have
changed since the VOTCA paper in 2009, which are:
-extrapolation function (quadratic -> exponential)
-interpolation function (cubic -> akima spline)
-shift (none -> minimun to 0)
Cheers,
Christoph
> --
> You received this message because you are subscribed to the Google Groups "votca" group.
> To post to this group, send email to vo...@googlegroups.com.
> To unsubscribe from this group, send email to votca+un...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/votca?hl=en.
>
--
Christoph Junghans
Web: http://www.compphys.de
Tinashe
I have noted the following error messages:
(Error 1)
$ awk -v kbt=1.6629 '{print $1, -kbt*log($2/$1/$1)}'
bond.dist.resample > bond.pot.ib
awk: (FILENAME=bond.dist.resample FNR=1) fatal: division by zero
attempted
Using a file (bond.pot.ib) which I obtained from csg_boltzmann:
(Error 2)
$ csg_call --options bond.xml --ia-type bond convert_potential gromacs
bond.pot.ib table_b1.xvg
Running subscript 'potential_to_gromacs.sh bond.pot.ib
table_b1.xvg'(from tags convert_potential gromacs)
Convert bond.pot.ib to table_b1.xvg
Running critical command 'mktemp bond.pot.smooth.XXXXX'
Running critical command 'csg_resample --in bond.pot.ib --out
bond.pot.smooth.18084 --grid 0:0.002: --comment Created on Wed Feb 15
17:43:06 CET 2012
called from potential_to_gromacs.sh, version 1.2.1 hgid: 88d17d52748a
settings file: /tmp/tinashe/votca-tutorials-1.2.1/test_propane/
atomistic_propane/bond.xml
working directory: /tmp/tinashe/votca-tutorials-1.2.1/test_propane/
atomistic_propane'
wrong range format, use min:step:max
###################################################################################################################################################################################################
#
ERROR:
#
# critical: 'csg_resample --in bond.pot.ib --out bond.pot.smooth.18084
--grid 0:0.002:
called from potential_to_gromacs.sh, version 1.2.1 hgid: 88d17d52748a
settings file: /tmp/tinashe/votca-tutorials-1.2.1/test_propane/
atomistic_propane/bond.xml
working directory: /tmp/tinashe/votca-tutorials-1.2.1/test_propane/
atomistic_propane' failed #
#
##################################################################################################################################################################################################Terminated
The following is the bond.xml file (I have also changed the table_end
value to 3.0 but still get the same error):
$ cat bond.xml
<inverse>
<program>gromacs</program>
<gromacs>
<pot_max>1e8</pot_max>
</gromacs>
</inverse>
</cg>
>
> --
> Christoph Junghans
> Web:http://www.compphys.de- Hide quoted text -
>
> - Show quoted text -
Christoph Junghans
> (Error 1)
> $ awk -v kbt=1.6629 '{print $1, -kbt*log($2/$1/$1)}'
> bond.dist.resample > bond.pot.ib
> awk: (FILENAME=bond.dist.resample FNR=1) fatal: division by zero
> attempted
>
> Using a file (bond.pot.ib) which I obtained from csg_boltzmann:
( $csg_resample --in bond.dist.new --out bond.dist.resample --grid
0.154:0.001:0.179)
As it crops away all point with y=0 and also x=0, which is why the
command fails in your case.
Strange. Which version of votca is that? There was a bug in that
script in version 1.2.0, but you seem to have 1.2.1.
Can you please post the output of
$csg_property --file bond.xml --path cg.inverse.gromacs.table_end
--short --print .
and
$csg_call --cat convert_potential gromacs
Christoph
> For more options, visit this group at http://groups.google.com/group/votca?hl=en.
Tinashe
When making the file "bond.dist.resample", the script automatically
puts a "#" on the first line and this is what was causing the error:
attempted
as output:
$ csg_property --file bond.xml --path cg.inverse.gromacs.table_end --
short --print .
0.3
Your second request is as follows:
$ csg_call --cat convert_potential gromacs
#!/bin/bash
#
# Copyright 2009-2011 The VOTCA Development Team (http://
www.votca.org)
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
if [[ $1 = "--help" ]]; then
cat <<EOF
${0##*/}, version %version%
This script is a wrapper to convert a potential to gromacs
Usage: ${0##*/} [input] [output]
EOF
exit 0
fi
if [[ -n $1 ]]; then
name="${1%%.*}"
input="$1"
shift
else
name=$(csg_get_interaction_property name)
input="${name}.pot.cur"
fi
[[ -f $input ]] || die "${0##*/}: Could not find input file '$input'"
if [[ -n $1 ]]; then
output="$1"
shift
else
output="$(csg_get_interaction_property inverse.gromacs.table)"
fi
echo "Convert $input to $output"
zero=0
tabtype="$(csg_get_interaction_property bondtype)"
#do this with --allow-empty to avoid stoping if calling from csg_call
[[ $(csg_get_property --allow-empty cg.inverse.method) = "tf" ]] &&
tabtype="thermforce"
if [[ $tabtype = "non-bonded" || $tabtype = "C6" || $tabtype =
"C12" ]]; then
tablend="$(csg_get_property --allow-empty
cg.inverse.gromacs.table_end)"
mdp="$(csg_get_property cg.inverse.gromacs.mdp "grompp.mdp")"
if [[ -f ${mdp} ]]; then
echo "Found setting file '$mdp' now trying to check options in
there"
rlist=$(get_simulation_setting rlist)
tabext=$(get_simulation_setting table-extension)
# if we have all 3 numbers do this checks
tabl=$(csg_calc "$rlist" + "$tabext")
[[ -n $tablend ]] && csg_calc "$tablend" "<" "$tabl" && \
die "${0##*/}: Error table is shorter then what mdp file ($mdp)
needs, increase cg.inverse.gromacs.table_end in setting file.\nrlist
($rlist) + tabext ($tabext) > cg.inverse.gromacs.table_end ($tablend)"
[[ -z $tablend ]] && tablend=$(csg_calc "$rlist" + "$tabext")
elif [[ -z $tablend ]]; then
die "${0##*/}: cg.inverse.gromacs.table_end was not defined in xml
seeting file"
fi
elif [[ $tabtype = "bonded" || $tabtype = "thermforce" ]]; then
tablend="$(csg_get_property cg.inverse.gromacs.table_end)"
elif [[ $tabtype = "angle" ]]; then
tablend=180
elif [[ $tabtype = "dihedral" ]]; then
zero="-180"
tablend=180
fi
gromacs_bins="$(csg_get_property cg.inverse.gromacs.table_bins)"
comment="$(get_table_comment $input)"
smooth="$(critical mktemp ${name}.pot.smooth.XXXXX)"
critical csg_resample --in ${input} --out "$smooth" --grid "${zero}:$
{gromacs_bins}:${tablend}" --comment "$comment"
extrapol="$(critical mktemp ${name}.pot.extrapol.XXXXX)"
tshift="$(critical mktemp ${name}.pot.shift.XXXXX)"
if [[ $tabtype = "non-bonded" || $tabtype = "C6" || $tabtype =
"C12" ]]; then
extrapol1="$(critical mktemp ${name}.pot.extrapol2.XXXXX)"
do_external table extrapolate --function exponential --avgpoints 5 --
region left "${smooth}" "${extrapol1}"
do_external table extrapolate --function constant --avgpoints 1 --
region right "${extrapol1}" "${extrapol}"
do_external pot shift_nonbonded "${extrapol}" "${tshift}"
elif [[ $tabtype = "thermforce" ]]; then
do_external table extrapolate --function constant --avgpoints 5 --
region leftright "${smooth}" "${extrapol}"
do_external pot shift_bonded "${extrapol}" "${tshift}"
else
do_external table extrapolate --function exponential --avgpoints 5 --
region leftright "${smooth}" "${extrapol}"
do_external pot shift_bonded "${extrapol}" "${tshift}"
fi
potmax="$(csg_get_property --allow-empty cg.inverse.gromacs.pot_max)"
[[ -n ${potmax} ]] && potmax="--max ${potmax}"
do_external convert_potential xvg ${potmax} --type "${tabtype}" "$
{tshift}" "${output}"
To create the angle and dihedral potentials, are the following
commands correct (I can email you the plots to have a more complete
view)?:
sed -e 's/$/ i/' BAB-angle.pot.ib | tac | sed -e '1,2d' | tac >
BAB.cut
csg_call table smooth BAB.cut BAB.smooth
csg_resample --in BAB.smooth --out BAB.refined --grid
0::0.001:3.141592654
csg_call table extrapolate --function sasha --region left BAB.refined
BAB.refined
csg_call table extrapolate --function sasha --region right BAB.refined
BAB.refined
awk '{print $1/3.141592654*180.0,$2}' BAB.refined > BAB.pot.cur
csg_call --ia-type angle --ia-name BAB --options convert.xml
convert_potential gromacs
sed -e '1d' -e '$d' BBAB-dihedral.pot.ib > BBAB.cut
#csg_call table smooth BBAB.cut BBAB.smooth <--- this command fails
and I renamed BBAB.cut to BBAB.smooth to proceed to the next command
csg_resample --in BBAB.smooth --out BBAB.refined --grid
-3.141592654::0.001:3.141592654
csg_call table extrapolate --function sasha --region left BBAB.refined
BBAB.refined
csg_call table extrapolate --function sasha --region right
BBAB.refined BBAB.refined
awk '{print $1/3.141592654*180.0,$2}' BBAB.refined > BBAB.pot.cur
csg_call --ia-type dihedral --ia-name BBAB --options convert.xml
convert_potential gromacs
Thanks.
Tinashe
---
On Feb 15, 6:37 pm, Christoph Junghans <jungh...@votca.org> wrote:
> Am 15. Februar 2012 09:58 schrieb Tinashe <tinashe.nd...@gmail.com>:> (Error 1)
> > ERROR:
> > #
> > # critical: 'csg_resample --in bond.pot.ib --out bond.pot.smooth.18084
> > --grid 0:0.002:
> > called from potential_to_gromacs.sh, version 1.2.1 hgid: 88d17d52748a
> > settings file: /tmp/tinashe/votca-tutorials-1.2.1/test_propane/
> > atomistic_propane/bond.xml
> > working directory: /tmp/tinashe/votca-tutorials-1.2.1/test_propane/
> > atomistic_propane' failed #
> > #
> >> - Show quoted text -
>
> > --
> > You received this message because you are subscribed to the Google Groups "votca" group.
> > To post to this group, send email to vo...@googlegroups.com.
> > To unsubscribe from this group, send email to votca+un...@googlegroups.com.
> > For more options, visit this group athttp://groups.google.com/group/votca?hl=en.
>
> --
> Christoph Junghans
> Web:http://www.compphys.de- Hide quoted text -
>
Christoph Junghans
> Hi Christoph,
>
> When making the file "bond.dist.resample", the script automatically
> puts a "#" on the first line and this is what was causing the error:
> awk: (FILENAME=bond.dist.resample FNR=1) fatal: division by zero
> attempted
away line starting with #
$ awk -v kbt=1.6629 '/^[^#]/{print $1, -kbt*log($2/$1/$1)}' >
bond.dist.resample > bond.pot.ib
/^[^#]/ means lines starting with any character but #.
> The following command gives the tabel_end value specified in bond.xml
> as output:
> $ csg_property --file bond.xml --path cg.inverse.gromacs.table_end --
> short --print .
> 0.3
>
> Your second request is as follows:
> $ csg_call --cat convert_potential gromacs
> #!/bin/bash
> ...
Both things look ok, which makes it hard to tell where this error comes from.
I would recommend to update to VOTCA 1.2.2 to see if the problem persists.
>
> To create the angle and dihedral potentials, are the following
> commands correct (I can email you the plots to have a more complete
> view)?:
Looks good!
>
> sed -e 's/$/ i/' BAB-angle.pot.ib | tac | sed -e '1,2d' | tac >
> BAB.cut
> csg_call table smooth BAB.cut BAB.smooth
> csg_resample --in BAB.smooth --out BAB.refined --grid
> 0::0.001:3.141592654
> csg_call table extrapolate --function sasha --region left BAB.refined
> BAB.refined
> csg_call table extrapolate --function sasha --region right BAB.refined
> BAB.refined
> awk '{print $1/3.141592654*180.0,$2}' BAB.refined > BAB.pot.cur
> csg_call --ia-type angle --ia-name BAB --options convert.xml
> convert_potential gromacs
>
> sed -e '1d' -e '$d' BBAB-dihedral.pot.ib > BBAB.cut
> #csg_call table smooth BBAB.cut BBAB.smooth <--- this command fails
You don't necessarily have to smooth the potential, only do that if is rough.
What was the error message here? (Update to Votca 1.2.2 first)
Christoph
> For more options, visit this group at http://groups.google.com/group/votca?hl=en.