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