Repository navigation
Expand file tree
/
Copy pathbar
More file actions
executable file
·99 lines (89 loc) · 2.91 KB
/
Copy pathbar
File metadata and controls
executable file
·99 lines (89 loc) · 2.91 KB
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
#!/usr/bin/perl
# BAR is a Perl script for analyzing the cumulative Helmholtz free
# energy change of a forward and backward whole-step perturbation
# from a MD trajectory using the Bennet Acceptance Ratio method
#
# The script will look for two key files and two parameter files
# associated with each simulation. Requires two sets of trajectory
# snapshots in different folders, along with separate parameter
# and key files in each folder.
#
# Written by Chris King, Ren Lab, University of Texas at Austin
#command-line arguments
if(@ARGV){
$struct=$ARGV[0];
$firstframe=$ARGV[1];
$lastframe=$ARGV[2];
@keys=($ARGV[3],$ARGV[4]);
$C=$ARGV[5];
@folders=($ARGV[6],$ARGV[7]);
}
#or get arguments from prompts
else{
print "Input Name of *.xyz Structure (No Extension!): ";
chomp($struct=<STDIN>);
print "First MD Frame Number (No Leading Zeros!): ";
chomp($firstframe=<STDIN>);
print "Last MD Frame Number (No Leading Zeros!): ";
chomp($lastframe=<STDIN>);
print "Input Two .key Filenames (with extensions) Separated by Spaces in Increasing Order of Lambda: ";
@keys= split ' ', <STDIN>;
print "Input starting energy: ";
chomp($C=<STDIN>);
print "Input Two Folder Names Separated by Spaces in Increasing Order of Lambda: ";
@folders= split ' ', <STDIN>;
}
#gas constant in kcal/mol and temperature
$R=.0019872066;
$T=298.0;
#compute avg exp terms, uses analysis subroutine
$RT=$R*$T;
$framenumber=$firstframe;
while($framenumber<=$lastframe){
if($framenumber<10){
$frame="$struct"."."."00"."$framenumber";}
elsif($framenumber<100){
$frame="$struct"."."."0"."$framenumber";}
else {$frame="$struct"."."."$framenumber";}
$u10[$framenumber]= &nrg($folders[1],$frame,$keys[0]);
$u11[$framenumber]= &nrg($folders[1],$frame,$keys[1]);
$u00[$framenumber]= &nrg($folders[0],$frame,$keys[0]);
$u01[$framenumber]= &nrg($folders[0],$frame,$keys[1]);
$framenumber+=1;
}
$framenumber=$firstframe;
$old=$C;
print "$old\n";
$new=&loop($old);
print "$new\n";
$diff=abs($new-$old);
while($diff >= 0.01) {
$old=$new;
$new=&loop($old);
$diff=abs($new-$old);
$bailout+=1;
if($bailout==100){
print "The Simulation Completed 100 Iterations without Converging.\n";
system "time";
}
print "$new\n";
}
sub loop {
$framenumber=$firstframe;
$num=0;
$denom=0;
while($framenumber<=$lastframe){
$num+=1/(1+exp(($u10[$framenumber]-$u11[$framenumber]+$_[0])/$RT));
$denom+=1/(1+exp(($u01[$framenumber]-$u00[$framenumber]-$_[0])/$RT));
$framenumber+=1;
}
$dA=$RT*log($num/$denom)+$_[0];
return $dA;
}
#free energy data extraction subroutine invokes tinker/analyze.x
sub nrg{
chdir "$_[0]";
my @analyze_string = split ' ', `analyze.x $_[1] E -k $_[2] | grep 'Total'`;
chdir "\.\.";
return $analyze_string[4];
}