Perl Script: Amin
This is a simple perl script
that will take a AMBER energy minimisation request issued by
Swiss-PdbViewer and process it. It requires that you have previously
modified your .cshrc login script to include this
In any case, do not forget to alter the paths (displayed in red) to point at the proper location on your
system.
Note: this script will use the 1991 force field parameters, and can be easily modified to run with the 94 parameters.
Note: This script also requires that you provide a default file for edit and parm. They must be in the directory from where you launch your minimisation job. If none are present, default files are automatically copied from a directory $amber/default. As this directory is not present in the standard distribution, you must create it and put a default file in. Default files are provided below.
|
Default "edtin" file (7 lines) to be put in directory $amber/default/ |
default edtin script
0 0 0 0
XYZ
OMIT
XRAY
0 0 0 0 0
QUIT
|
|
Default "prmin" file (5 lines) to be put in directory $amber/default/ |
default prmin
BIN FOR STDA
0
|
#!/usr/central/bin/perl # split PDB file from AMBER command files # and preform the requested Energy Minimization # # version 1.0 # NG/97 $amber = $ENV{"amber"}; # ------------------------------------- # -------- retrieve arguments --------- # ------------------------------------- if ($#ARGV == -1) { print "v1.0 usage: Amin [-r] [-p] filename\n"; exit; } while ($#ARGV > 0) { if ($ARGV[0] eq "-p") { $parallel++; shift(@ARGV); } if ($ARGV[0] eq "-r") { $cleanup++; shift(@ARGV); } } $pdbfile = @ARGV[0]; if ($parallel) { $sander = $ENV{"pamber"}."/sander"; } else { $sander = "sander"; } # ----------------------------------------- # -------- prepare some filenames --------- # ----------------------------------------- ($minifile) = split(/\./,$pdbfile); $infile = $$.".pdb"; $protonatedfile = $$.".H.pdb"; $linkfile = $$.".lnkin"; $sanderfile = $$.".sanderin"; $sanderoutfile = $$.".sanderout"; $sandermdinfo = $$.".mdinfo"; $tmpfile = $$.".tmp"; $minimizedXYZ = $$.".xyz"; $minimizedPDB = $$. "._E.pdb"; $minimizedfile = $minifile . "_E.pdb"; # --------------------------------------------------------------- # ------ Create an input file with proper line terminators ------ # ------ (for mac translate CR in LF; for pc remove CR) ------ # --------------------------------------------------------------- # ---- figure out which kind of file we deal with --- open (IN,"<$pdbfile") or die "** ERROR: file $pdbfile doesn't exist\n"; (@line) = split '', <IN>; $nbchar = $#line; if ($nbchar > 200) { $nbchar = 200; } for ($i=0; $i <= $nbchar; $i++) { if ($line[$i] eq "\n") { $lf = 1; } elsif ($line[$i] eq "\r") { $cr = 1; } } close IN; # ----- do the appropriate conversion --- open (IN,"<$pdbfile") or die; open (TMP,">$tmpfile") or die; while(<IN>) { if ($cr && $lf) { ~ s/[\r]//; } else { ~ tr/[\r]/[\n]/; } print TMP $_; } close TMP; close IN; if (-e "%$pdbfile") { unlink "%$pdbfile"; } # ----------------------------------------------------------------------- # ------- Split input file in PDB, LINK and SANDER (very crude) --------- # ----------------------------------------------------------------------- open (TMP,"<$tmpfile") or die; open (PDB,">$infile") or die; open (LINK,">$linkfile") or die; open (SANDER,">$sanderfile") or die; while(<TMP>) { print PDB $_; if (substr($_,0,3) eq "END") { while(<TMP>) { if (substr($_,0,4) eq "*EOF") { (@CHAINS) = split '', <TMP>; while(<TMP>) { print SANDER $_; } } else{ print LINK $_; } } } } close SANDER; close LINK; close PDB; close TMP; # ------------------------------------------------- # --------- Minimize (at least try to...) --------- # ------------------------------------------------- unless (-e "edtin") { system("cp $amber/default/edtin ."); print "edtin file not found, taking default\n" } unless (-e "prmin") { system("cp $amber/default/prmin ."); print "prmin file not found, taking default\n" } print "** PREPARING **\n"; system("protonate -d $amber/dat/PROTON_INFO < $infile > $protonatedfile"); system("link -i $linkfile -l $$.lnkbin -o $$.lnkout -p $amber/dat/db4.dat"); system("edit -i edtin -l $$.lnkbin -e $$.edtbin -o $$.edtout -pi $protonatedfile"); system("parm -i prmin -o $$.prmout -e $$.edtbin -f $amber/dat/parm91.dat -c $$.prmcrd -p $$.prmtop"); print "** MINIMIZING **\n"; if ($parallel) { print "(using parallel version)\n"; } system("$sander -i $sanderfile -o $sanderoutfile -c $$.prmcrd -p $$.prmtop -r $minimizedXYZ -inf $sandermdinfo"); if (-e $minimizedXYZ) { # WARNING: ambpdb must take file named "prmtop" as input. The way I deal with this # can cause problems if two processes are running in parallel system("mv $$.prmtop prmtop"); system("ambpdb < $minimizedXYZ > $minimizedPDB"); system("mv prmtop $$.prmtop "); # --- Add energy minimisation results as comments at top of the pdb file open(ENERGY, "< $sandermdinfo") or die; open (RESULT,">$minimizedfile") or die; print RESULT "// Structure minimized with AMBER\n"; while(<ENERGY>) { print RESULT "//".$_; } close ENERGY; # --- Write energy minimized file and add CHAIN label $i = 0; open(MINIMIZED, "< $minimizedPDB") or die; while(<MINIMIZED>) { if (substr($_,0,4) eq "ATOM") { print RESULT substr($_,0,21).$CHAINS[$i].substr($_,22); } elsif (substr($_,0,3) eq "TER") { print RESULT $_; $i++; } } close MINIMIZED; close RESULT; print "** DONE **\n"; if ($cleanup) { system("\\rm $$.*"); } } else { print "** ERROR ** (process id = ".$$.")\n"; } exit;