Hi Christoph, thx again for your support this great tool. Im glad I didnt give up. Meanwhile I've tested my new cg-potential in a tensile test of 40 chains of PP each having 2400 repeats units. The result is surprisignly close to the experiment (simulation is still running) and i can use a timestep of 30 fs.
I surley can share my little code-modifications. But its very sloppy and I just built a little "workaround". However, I've appended the modified files and if you replace the original scripts, votca will only work with lammps and no pressure correction is implemented yet. So do a backup of those files. csg_stat is replaced by a simple bash script. The rdfs and distributions are calculated within lammps as defined in the lammps-input-script. The lammps-script calls another bash-script which converts the tables from gromacs-units to units-real (see lammps doc:distance=Angstrom, Energy=kcal/mole). This file has to be added into the <filelist> -section in your settings.xml. The options for the histrograms need to be adjusted in the lammps-input script. The values defined in settings.xml are ignored. Keep in mind: Im not a professional, this is only a little hack. Its not a full implementation and im not a about to learn perl to do so. Its working for me and maybe it will work for somebody else too.. regards, frank. -- 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 [email protected]. To post to this group, send email to [email protected]. Visit this group at http://groups.google.com/group/votca. For more options, visit https://groups.google.com/d/optout.
csg_stat
Description: Binary data
#! /usr/bin/perl -w # # Copyright 2009-2014 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. # use strict; ( my $progname = $0 ) =~ s#^.*/##; my $usage="Usage: $progname [OPTIONS] <in> <derivatives_in> <out>"; my $type="non-bonded"; my $sim_prog="none"; while ((defined ($ARGV[0])) and ($ARGV[0] =~ /^-./)) { if (($ARGV[0] !~ /^--/) and (length($ARGV[0])>2)){ $_=shift(@ARGV); #short opt having agruments examples fo if ( $_ =~ /^-[fo]/ ) { unshift(@ARGV,substr($_,0,2),substr($_,2)); } else{ unshift(@ARGV,substr($_,0,2),"-".substr($_,2)); } } if (($ARGV[0] eq "-h") or ($ARGV[0] eq "--help")){ print <<EOF; $progname, version %version% This script converts csg potential files to the tab format (as read by espresso or lammps or dlpoly). In addition, it does some magic tricks: - shift the potential, so that it is zero at the cutoff $usage Allowed options: -h, --help show this help message --type XXX change the type of xvg table Default: $type --header XXX Write a special simulation programm header Examples: * $progname --type non-bonded table.in table_b0.xvg EOF exit 0; } elsif ($ARGV[0] eq "--type"){ shift(@ARGV); $type = shift(@ARGV); } elsif ($ARGV[0] eq "--header"){ shift(@ARGV); $sim_prog = shift(@ARGV); } else{ die "Unknown option '".$ARGV[0]."' !\n"; } } die "$progname: conversion of ${type} interaction to generic tables is only implemented for dlpoly and lammps!\n" unless (($type eq "non-bonded")||($sim_prog eq "dlpoly")||($sim_prog eq "lammps")); die "3 parameters are necessary\n" if ($#ARGV<2); use CsgFunctions; my $in_pot="$ARGV[0]"; my $in_deriv_pot="$ARGV[1]"; my $outfile="$ARGV[2]"; my @r; my @r_repeat; my @pot; my @pot_deriv; my @flag; my @flag_repeat; #cutoff is last point (readin_table($in_pot,@r,@pot,@flag)) || die "$progname: error at readin_table\n"; (readin_table($in_deriv_pot,@r_repeat,@pot_deriv,@flag_repeat)) || die "$progname: error at readin_table\n"; if ($type eq "non-bonded"){ #shift potential so that it is zero at cutoff for (my $i=0;$i<=$#r;$i++){ $pot[$i]-=$pot[$#r]; } } #hopefully shift for bonded interactions have been done outside open(OUTFILE,"> $outfile") or die "saveto_table: could not open $outfile\n"; if ($sim_prog eq "espresso") { # espresso specific header - no other starting comments printf(OUTFILE "#%d %f %f\n", $#r+1, $r[0],$r[$#r]); for(my $i=0;$i<=$#r;$i++){ printf(OUTFILE "%15.10e %15.10e %15.10e\n",$r[$i],($r[$i]>0)?-$pot_deriv[$i]/$r[$i]:-$pot_deriv[$i], $pot[$i]); } } elsif ($sim_prog eq "lammps") { #achtung: ich konervierte alles in units real einheiten if ($type eq "non-bonded"){ printf(OUTFILE "VOTCA\n"); printf(OUTFILE "N %i R %f %f\n\n",$#r+1,$r[0],$r[$#r]); #printf(OUTFILE "N %i R %f %f\n\n",$#r+1,$r[0]*10,$r[$#r]*1r 0); for(my $i=0;$i<=$#r;$i++){ #printf(OUTFILE "%i %15.10e %15.10e %15.10e\n",$i+1,$r[$i]*10, $pot[$i]*0.239, -$pot_deriv[$i]*0.0239); printf(OUTFILE "%i %15.10e %15.10e %15.10e\n",$i+1,$r[$i],$pot[$i],-$pot_deriv[$i]); } } elsif ( $type eq "bond" ) { printf(OUTFILE "VOTCA\n"); printf(OUTFILE "N %i\n\n",$#r+1); for(my $i=0;$i<=$#r;$i++){ #nm -> Angs: $r[$i]*10.0 printf(OUTFILE "%i %12.5e %15.7e %15.7e\n",$i+1,$r[$i],$pot[$i],-$pot_deriv[$i]*$r[$i]); } } elsif ( $type eq "angle" || $type eq "dihedral" ) { printf(OUTFILE "VOTCA\n"); printf(OUTFILE "N %i\n\n",$#r+1); my $RadToDegree=180/3.14159265359; for(my $i=0;$i<=$#r;$i++){ #rad -> degree: $r[$i]*$RadToDegree, and $pot_deriv[$i]/$RadToDegree printf(OUTFILE "%i %12.5e %15.7e %15.7e\n",$i+1,$r[$i], $pot[$i], -$pot_deriv[$i]/$RadToDegree); #printf(OUTFILE "%i %12.5e %15.7e %15.7e\n",$i+1,$r[$i]*$RadToDegree, $pot[$i]*0.239, -$pot_deriv[$i]*0.0239/$RadToDegree); } } } elsif ($sim_prog eq "dlpoly") { if ($type eq "non-bonded"){ # see dlpoly manual ngrid = cut/delta + 4 = $#r + 4 as table starts with delta (not 0) # number of lines int((ngrid+3)/4) for(my $i=0;$i<4*int(($#r+7)/4);$i++){ printf(OUTFILE "%15.7e",($i>$#r)?0:$pot[$i]); printf(OUTFILE "%s",($i%4==3)?"\n":" "); } for(my $i=0;$i<4*int(($#r+7)/4);$i++){ # no scaling factor needed 1 kJ/nm *nm = 1 (kJ/Angs)*Angs printf(OUTFILE "%15.7e",($i>$#r)?0:-$pot_deriv[$i]*$r[$i]); printf(OUTFILE "%s",($i%4==3)?"\n":" "); } } elsif ( $type eq "bond" ) { for(my $i=0;$i<=$#r;$i++){ #nm -> Angs: $r[$i]*10.0 printf(OUTFILE "%12.5e %15.7e %15.7e\n",$r[$i]*10.0, $pot[$i], -$pot_deriv[$i]*$r[$i]); } } elsif ( $type eq "angle" || $type eq "dihedral" ) { my $RadToDegree=180.0/3.14159265359; for(my $i=0;$i<=$#r;$i++){ #rad -> degree: $r[$i]*$RadToDegree, and $pot_deriv[$i]/$RadToDegree printf(OUTFILE "%12.5e %15.7e %15.7e\n",$r[$i]*$RadToDegree, $pot[$i], -$pot_deriv[$i]/$RadToDegree); } } else { #should never happen die "$progname: tabulated potentials/forces for dlpoly $type not implemented\n"; } } elsif ($sim_prog eq "gromacs") { printf(OUTFILE "#This is just a failback, for using different columns use table_to_xvg.pl instead!\n"); for(my $i=0;$i<=$#r;$i++){ printf(OUTFILE "%15.10e %15.10e %15.10e %15.10e %15.10e %15.10e %15.10e\n",$r[$i], ,0,0,0,0,$pot[$i], -$pot_deriv[$i]); } } else { #generic for espressopp / hoomd-blue for(my $i=0;$i<=$#r;$i++){ printf(OUTFILE "%15.10e %15.10e %15.10e\n",$r[$i], $pot[$i], -$pot_deriv[$i]); } } close(OUTFILE) or die "Error at closing $outfile\n";
#!/bin/bash
cp A-A.pot tmp.A-A.pot
awk '{
if(NR==2){$4=$4*10;$5=$5*10;}
if(NF==4){$2=$2*10;$3=$3*0.239;$4=$4*0.0239*0.239;}
print $0;
}' tmp.A-A.pot > A-A.pot.lmps
cp bond.pot tmp.bond.pot
awk '{
if(NF==4){$2=$2*10;$3=$3*0.239;$4=$4*0.0239*0.239;}
print $0;
}' tmp.bond.pot > bond.pot.lmps
cp angle.pot tmp.angle.pot
awk 'BEGIN{rad2deg=180/3.14159265358979323;}{
if(NF==4){$2=$2;$3=$3*0.239;$4=$4*0.0239*0.239/rad2deg;}
print $0;
}' tmp.angle.pot > angle.pot.lmps
my.in
Description: Binary data
