1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
|
#! /usr/bin/env -S perl -w
#
# Copyright 2009-2017 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#^.*/##;
if (defined($ARGV[0])&&("$ARGV[0]" eq "--help")){
print <<EOF;
$progname, version %version%
This script calcs dU out of two rdfs with the rules of inverse boltzmann
In addition, it does some magic tricks:
- do not update if one of the two rdf is undefined
Usage: $progname target_rdf new_rdf cur_pot outfile kBT
EOF
exit 0;
}
die "5 parameters are necessary\n" if ($#ARGV<4);
use CsgFunctions;
my $pref="$ARGV[4]";
my $aim_rdf_file="$ARGV[0]";
my @r_aim;
my @rdf_aim;
my @flags_aim;
(readin_table($aim_rdf_file,@r_aim,@rdf_aim,@flags_aim)) || die "$progname: error at readin_table\n";
my $cur_rdf_file="$ARGV[1]";
my @r_cur;
my @rdf_cur;
my @flags_cur;
(readin_table($cur_rdf_file,@r_cur,@rdf_cur,@flags_cur)) || die "$progname: error at readin_table\n";
my $cur_pot_file="$ARGV[2]";
my @pot_r_cur;
my @pot_cur;
my @pot_flags_cur;
(readin_table($cur_pot_file,@pot_r_cur,@pot_cur,@pot_flags_cur)) || die "$progname: error at readin_table\n";
#should never happen due to resample, but better check
die "Different grids \n" if (($r_aim[1]-$r_aim[0]-$r_cur[1]+$r_cur[0]) > 0.0001);
die "Different start potential point \n" if (($r_aim[0]-$r_cur[0]) > 0.0001);
die "Different end potential point \n" if ( $#r_aim != $#r_cur );
my $outfile="$ARGV[3]";
my @dpot;
my @flag;
# return the index of largest value in an array
sub maxindex {
my $index = 0;
my $max = 0.0;
for my $i (0 .. $#_) {
if ($_[$i] > $max) {
$max = $_[$i];
$index = $i;
}
}
return $index;
}
my $ndx_max_rdf_cur = maxindex(@rdf_cur);
# go from ndx_max_rdf to end of table and fill the end with the last valid value encountered
my $value=0.0;
for (my $i=$ndx_max_rdf_cur;$i<=$#r_aim;$i++){
if (($rdf_aim[$i] > 1e-10) && ($rdf_cur[$i] > 1e-10)) {
$dpot[$i]=log($rdf_cur[$i]/$rdf_aim[$i])*$pref;
$flag[$i]="i";
} else {
$dpot[$i]=$value;
$flag[$i]="o";
}
if($pot_flags_cur[$i] =~ /[u]/) {
$dpot[$i]=$value;
$flag[$i]="o";
}
else {
$value=$dpot[$i];
}
}
# go from ndx_max_rdf to beginning of table and fill the end with the last valid value encountered
$value=0.0;
for (my $i=$ndx_max_rdf_cur-1;$i>=0;$i--){
if (($rdf_aim[$i] > 1e-10) && ($rdf_cur[$i] > 1e-10)) {
$dpot[$i]=log($rdf_cur[$i]/$rdf_aim[$i])*$pref;
$flag[$i]="i";
} else {
$dpot[$i]=$value;
$flag[$i]="o";
}
if($pot_flags_cur[$i] =~ /[u]/) {
$dpot[$i]=$value;
$flag[$i]="o";
}
else {
$value=$dpot[$i];
}
}
my $comment="#progname: aim_rdf=$aim_rdf_file cur_rdf=$cur_rdf_file cur_pot=$cur_rdf_file\n";
saveto_table($outfile,@r_aim,@dpot,@flag,$comment) || die "$progname: error at save table\n";
|
