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
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 # # 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 <[email protected]> wrote: > Am 15. Februar 2012 09:58 schrieb Tinashe <[email protected]>:> (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: > > Have you done the resampling like I did above. > ( $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. > > > > > > > > > (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 > > 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 > > > > > > > > > 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 > > <cg> > > <inverse> > > <program>gromacs</program> > > <gromacs> > > <pot_max>1e8</pot_max> > > <table_end>0.3</table_end> > > <table_bins>0.002</table_bins> > > </gromacs> > > </inverse> > > </cg> > > > On Feb 14, 6:49 pm, Christoph Junghans <[email protected]> wrote: > >> Hi Tinashe, > > >> 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 > > >> Am 13. Februar 2012 04:23 schrieb Tinashe <[email protected]>: > > >> > I have been trying to reproduce the table_*.xvg files in the propane > >> > 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. > > >> > -- > >> > You received this message because you are subscribed to the Google > >> > Groups "votca" group. > >> > To post to this group, send email to [email protected]. > >> > To unsubscribe from this group, send email to > >> > [email protected]. > >> > For more options, visit this group > >> > athttp://groups.google.com/group/votca?hl=en. > > >> -- > >> Christoph Junghans > >> Web:http://www.compphys.de-Hide quoted text - > > >> - 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 [email protected]. > > To unsubscribe from this group, send email to > > [email protected]. > > For more options, visit this group > > athttp://groups.google.com/group/votca?hl=en. > > -- > Christoph Junghans > Web:http://www.compphys.de- Hide quoted text - > > - Show quoted text -- Hide quoted text - > > - 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 [email protected]. To unsubscribe from this group, send email to [email protected]. For more options, visit this group at http://groups.google.com/group/votca?hl=en.
