1#!/usr/local/bin/perl -w 2# 3# Copyright 2009-2011 The VOTCA Development Team (http://www.votca.org) 4# 5# Licensed under the Apache License, Version 2.0 (the "License"); 6# you may not use this file except in compliance with the License. 7# You may obtain a copy of the License at 8# 9# http://www.apache.org/licenses/LICENSE-2.0 10# 11# Unless required by applicable law or agreed to in writing, software 12# distributed under the License is distributed on an "AS IS" BASIS, 13# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14# See the License for the specific language governing permissions and 15# limitations under the License. 16# 17 18use strict; 19 20( my $progname = $0 ) =~ s#^.*/##; 21my $usage="Usage: $progname [OPTIONS] <in> <out>"; 22my $type="non-bonded"; 23my $gmx_max=undef; 24 25while ((defined ($ARGV[0])) and ($ARGV[0] =~ /^-./)) 26{ 27 if (($ARGV[0] !~ /^--/) and (length($ARGV[0])>2)){ 28 $_=shift(@ARGV); 29 #short opt having agruments examples fo 30 if ( $_ =~ /^-[fo]/ ) { 31 unshift(@ARGV,substr($_,0,2),substr($_,2)); 32 } 33 else{ 34 unshift(@ARGV,substr($_,0,2),"-".substr($_,2)); 35 } 36 } 37 if (($ARGV[0] eq "-h") or ($ARGV[0] eq "--help")){ 38 print <<EOF; 39$progname, version %version% 40 41This script converts csg potential files to the xvg format. 42 43$usage 44 45Allowed options: 46-h, --help show this help message 47--type XXX change the type of xvg table 48 Default: $type 49--max MAX Replace all pot value bigger MAX by MAX 50 51 52Possible types: non-bonded (=C12), bond, C12, C6, CB, angle, dihedral 53 54Examples: 55* $progname --type bond table.in table_b0.xvg 56EOF 57 exit 0; 58 } 59 elsif ($ARGV[0] eq "--type"){ 60 shift(@ARGV); 61 $type = shift(@ARGV); 62 } 63 elsif ($ARGV[0] eq "--max"){ 64 shift(@ARGV); 65 $gmx_max = shift(@ARGV); 66 } 67 else{ 68 die "Unknown option '".$ARGV[0]."' !\n"; 69 } 70} 71 72die "2 parameters are necessary\n" if ($#ARGV<1); 73 74use CsgFunctions; 75 76my $infile="$ARGV[0]"; 77my $outfile="$ARGV[1]"; 78 79my @r; 80my @pot; 81my @flag; 82my $comments; 83(readin_table($infile,@r,@pot,@flag,$comments)) || die "$progname: error at readin_table\n"; 84 85if (defined($gmx_max)) { 86 #gromacs does not like VERY big numbers 87 for (my $i=0;$i<=$#r;$i++) { 88 $pot[$i]=$gmx_max if $pot[$i]>$gmx_max; 89 $pot[$i]=-$gmx_max if $pot[$i]<-$gmx_max; 90 } 91} 92 93my @force; 94 95#calc force 96for (my $i=1;$i<$#r;$i++){ 97 $force[$i]=-($pot[$i+1]-$pot[$i-1])/($r[$i+1]-$r[$i-1]); 98} 99if ( "$type" eq "dihedral" ) { 100 $force[0]=-($pot[1]-$pot[$#r-1])/($r[1]-$r[0]+$r[$#r]-$r[$#r-1]); 101 $force[$#r]=$force[0]; 102} else { 103 $force[0]=0; 104 $force[$#r]=0.0; 105} 106 107open(OUTFILE,"> $outfile") or die "saveto_table: could not open $outfile\n"; 108 109my $fmt=undef; 110my $begin=0; 111my $end=undef; 112if (( "$type" eq "non-bonded" ) or ("$type" eq "C12" )) { 113 $fmt=sprintf("%%15.10e %15.10e %15.10e %15.10e %15.10e %%15.10e %%15.10e\n",0,0,0,0); 114} 115elsif ( "$type" eq "C6" ){ 116 $fmt=sprintf("%%15.10e %15.10e %15.10e %%15.10e %%15.10e %15.10e %15.10e\n",0,0,0,0); 117} 118elsif ( "$type" eq "CB" ){ 119 $fmt=sprintf("%%15.10e %%15.10e %%15.10e %15.10e %15.10e %15.10e %15.10e\n",0,0,0,0); 120} 121elsif ( "$type" eq "bond" ){ 122 $fmt="%15.10e %15.10e %15.10e\n"; 123} 124elsif ( "$type" eq "angle" ){ 125 $fmt="%15.10e %15.10e %15.10e\n"; 126 $end=180; 127} 128elsif ( "$type" eq "dihedral" ){ 129 $fmt="%15.10e %15.10e %15.10e\n"; 130 $begin=-180; 131 $end=180; 132} 133else{ 134 die "$progname: Unsupported type of interatction: $type -> go and implement it\n"; 135} 136 137die "$progname: table for type $type should begin with $begin, but I found $r[0]\n" if(abs($begin-$r[0]) > 1e-3); 138die "$progname: table for type $type should end with $end, but I found $r[$#r]\n" if(($end) and (abs($end-$r[$#r]) > 1e-3)); 139 140print OUTFILE "$comments" if (defined($comments)); 141for(my $i=0;$i<=$#r;$i++){ 142 printf(OUTFILE "$fmt",$r[$i],$pot[$i], $force[$i]); 143} 144close(OUTFILE) or die "Error at closing $outfile\n"; 145 146