1#!/usr/local/bin/perl -w
2#
3# Copyright 2009-2017 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] kbint target_kbint outfile kBT min:step:max int_start:int_end ramp_factor";
22
23while ((defined ($ARGV[0])) and ($ARGV[0] =~ /^-./))
24{
25  if (($ARGV[0] !~ /^--/) and (length($ARGV[0])>2)){
26    $_=shift(@ARGV);
27    #short opt having agruments examples fo
28    if ( $_ =~ /^-[fo]/ ) {
29      unshift(@ARGV,substr($_,0,2),substr($_,2));
30    } else{
31      unshift(@ARGV,substr($_,0,2),"-".substr($_,2));
32    }
33  }
34  if (($ARGV[0] eq "-h") or ($ARGV[0] eq "--help")){
35     print <<END;
36$progname, version %version%
37
38This script calculates Kirkwood-Buff correction as described in:
39P. Ganguly, D. Mukherji, C. Junghans, N. F. A. van der Vegt,
40Kirkwood-Buff coarse-grained force fields for aqueous solutions,
41J. Chem. Theo. Comp., 8, 1802 (2012), doi:10.1021/ct3000958
42
43$usage
44
45Allowed options:
46-h, --help            Show this help message
47END
48    exit 0;
49  }else{
50    die "Unknown option '".$ARGV[0]."' !\n";
51  }
52}
53
54die "7 parameters are necessary\n" if ($#ARGV<6);
55
56use CsgFunctions;
57
58my $kbt=$ARGV[3];
59my @irange=split(/:/,$ARGV[5]);
60defined($irange[1]) || die "Not enough number in irange $ARGV[5], got ".($#irange+1)." need 2\n";
61my $int_start=$irange[0];
62my $int_stop=$irange[1];
63my $ramp_factor=$ARGV[6];
64
65my @range=split(/:/,$ARGV[4]);
66defined($range[2]) || die "Not enough number in range $ARGV[4], got ".($#range+1)." need 3\n";
67my $r_min=$range[0];
68my $r_ramp=$range[2];
69my $delta_r=$range[1];
70
71my $aim_kbint_file="$ARGV[0]";
72my @r_aim;
73my @kbint_aim;
74my @flags_aim;
75(readin_table($aim_kbint_file,@r_aim,@kbint_aim,@flags_aim)) || die "$progname: error at readin_table\n";
76
77my $cur_kbint_file="$ARGV[1]";
78my @r_cur;
79my @kbint_cur;
80my @flags_cur;
81(readin_table($cur_kbint_file,@r_cur,@kbint_cur,@flags_cur)) || die "$progname: error at readin_table\n";
82
83#should never happen due to resample, but better check
84die "Different grids \n" if (($r_aim[1]-$r_aim[0]-$r_cur[1]+$r_cur[0])>0.0001);
85die "Different start potential point \n" if (($r_aim[0]-$r_cur[0]) > 0.0001);
86die "Different end potential point \n" if ( $#r_aim != $#r_cur );
87
88my $j=0;
89my $avg_int=0;
90for (my $i=0;$i<=$#r_aim;$i++){
91  if (($r_aim[$i]>=$int_start) && ($r_aim[$i]<=$int_stop)) {
92     $avg_int+=$kbint_cur[$i]-$kbint_aim[$i];
93     $j++;
94  }
95}
96$avg_int/=$j;
97
98my $comment="#$progname: avg_int($int_start:$int_stop)=$avg_int ramp_factor=$ramp_factor r_ramp=$r_ramp\n";
99my @dpot;
100my @flag;
101for (my $i=0;$i<=$#r_aim;$i++){
102  if ($r_aim[$i]> $r_ramp) {
103    $dpot[$i]=0; #beyond r_ramp correction is 0
104  } else {
105    $dpot[$i]=($avg_int*$ramp_factor*(1.0-($r_aim[$i]/$r_ramp)))*$kbt;
106  }
107  $flag[$i]="i";
108}
109
110my $outfile="$ARGV[2]";
111saveto_table($outfile,@r_aim,@dpot,@flag,$comment) || die "$progname: error at save table\n";
112