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