#!/usr/bin/perl -w

# parse command line arguments
if (scalar(@ARGV) < 3) {printHelpAndExit();}
## ARGUMENTS:
## parmfile (amber topology file) e.g. /home/icre85/icre85227/projects/CALIX/1BXXX_PNO/1BXXX_PNO_wat.top
## ambercoor (amber restart file) e.g. /home/icre85/icre85227/projects/CALIX/1BXXX_PNO/RESTRT/md_40.rst
## name of NAMD restart file (root only)
my $top = shift @ARGV;
my $rst = shift @ARGV;
my $name = shift @ARGV;

## quick check to ensure that files exist and are of the right type
#
if(-f $top && ($top =~ /\.top$/i || $top =~ /\.par$/i) ){
	print STDERR "First file seems ok\n";
}else{
	die "Are you sure that the first file is an AMBER topology file?";
}
if(-f $rst && ($rst =~ /\.rst$/i || $rst =~ /\.crd$/i) ){
	print STDERR "Second file seems ok\n";
}else{
	die "Are you sure that the second file is an AMBER restart file?";
}

# Assign name to velocities output files
my $ini = rindex($rst,"/") + 1;
my $end = rindex($rst,".");
my $length = $end - $ini;
my $root = substr($rst, $ini, $length);
my $vel = $root."_vel.rst";
my $vel_pdb = $root."_vel.pdb";

## READ RST FILE TO EXTRACT (Box size + timestep) AND CREATE velocities files
## 
my ($box, $step) = read_restart($rst,$vel);
$step = sprintf "%i",$step;

## EJECUTE ambpdb TO GET A PDB FILE WITH VELOCITIES
my $oo = `which ambpdb 2>&1`;
if ($oo =~ /^\//){
	my $output = `ambpdb -p $top < $vel > $vel_pdb`;
}else{
	print STDERR "WARNING - ambpdb not found.\n";
	print STDERR "\tMake sure it is in your PATH or use the following commad to obtain a PDB file with velocities:\n";
	print STDERR "\n\tambpdb -p $top < $vel > $vel_pdb\n\n";
}

## GUESS PME Grid size - should have only small integer factors (2, 3 and 5).
$pme = guess_pme($box);

## OUTPUT NAMD configuration file
my $conf = $name.".conf";
open(FC,">$conf") || die "Cannot open file $conf: $!";
printf FC "%s\n",q(# Definition of truncated octahedron box);
printf FC "%s\n",qq(set a $box ;      # aresta; taken from the Input CRD/RST file);
printf FC "%s\n",q(set p 109.471220634);
printf FC "%s\n",q(set pi 3.14159265358979323);
printf FC "%s\n",q(set cosp [expr cos($p * $pi/180)]);
printf FC "%s\n",q(set sinp [expr sin($p * $pi/180)]);
printf FC "%s\n",q(set ori [expr $a/2]);
printf FC "%s\n",q();
printf FC "%s\n",q(set Bx [expr $a*$cosp]);
printf FC "%s\n",q(set By [expr $a*$sinp]);
printf FC "%s\n",q(set Cx $Bx);
printf FC "%s\n",q(set Cy [expr $a*$cosp*(1-$cosp)/$sinp]);
printf FC "%s\n",q(set Cz [expr sqrt($a*$a - $Cx*$Cx - $Cy*$Cy)]);
printf FC "%s\n",q();
printf FC "%s\n",q(cellBasisVector1 $a 0 0);
printf FC "%s\n",q(cellBasisVector2 $Bx $By 0);
printf FC "%s\n",q(cellBasisVector3 $Cx $Cy $Cz);
printf FC "%s\n",q(cellOrigin $ori $ori $ori);
printf FC "%s\n",q(# END OF BOX DEFINITION);
printf FC "%s\n",q();
printf FC "%s\n",q(set temperature 300);
printf FC "%s\n",q(rigidBonds     all);
printf FC "%s\n",q();
printf FC "%s\n",q(pairlistdist   13.5);
printf FC "%s\n",q();
printf FC "%s\n",q(rigidTolerance 0.00001  # Default is  0.00001);
printf FC "%s\n",q(stepspercycle  10    # Num of timesteps between atom reassignment in non-bonded force evaluation);
printf FC "%s\n",q(timestep       2     # in unit of fs (This is default. Corresponds to dt in sander));
printf FC "%s\n",q(outputEnergies 1000  # Energy output frequency (corresponds to ntpr in sander));
printf FC "%s\n",q(restartfreq    500   # Restart file frequency (correspond to ntwr in sander));
printf FC "%s\n",q(DCDfreq        1000  # Trajectory file frequency (corresponds to ntwx in sander));
printf FC "%s\n",q();
printf FC "%s\n",q(# Non-bonded interaction parameter);
printf FC "%s\n",q(cutoff         9);
printf FC "%s\n",q(switching      off  # Not used when the Amber FF is used);
printf FC "%s\n",q();
printf FC "%s\n",q(# PME);
printf FC "%s\n",q(wrapAll on);
printf FC "%s\n",q(wrapNearest on);
printf FC "%s\n",q(PME on);
printf FC "%s\n",qq(PMEGridsizeX $pme);
printf FC "%s\n",qq(PMEGridsizeY $pme);
printf FC "%s\n",qq(PMEGridsizeZ $pme);
printf FC "%s\n",q();
printf FC "%s\n",q(# Pressure controlling);
printf FC "%s\n",q(useGroupPressure yes);
printf FC "%s\n",q(langevinPiston        on);
printf FC "%s\n",q(langevinPistonTarget  1.01325      ;# pressure in bar -> 1 atm);
printf FC "%s\n",q(langevinPistonPeriod  100.         ;# oscillation period around 100 fs);
printf FC "%s\n",q(langevinPistonDecay   50.          ;# oscillation decay time of 50 fs);
printf FC "%s\n",q(langevinPistonTemp    $temperature ;# coupled to heat bath);
printf FC "%s\n",q();
printf FC "%s\n",q(# Constant Temperature Control);
printf FC "%s\n",q(langevin            on            ;# langevin dynamics);
printf FC "%s\n",q(langevinDamping     5.            ;# damping coefficient of 5/ps);
printf FC "%s\n",q(langevinTemp        $temperature  ;# random noise at this level);
printf FC "%s\n",q(langevinHydrogen    no            ;# don't couple bath to hydrogens);
printf FC "%s\n",q();
printf FC "%s\n",q(# FF);
printf FC "%s\n",q(amber               on  # Specify this is AMBER force field);
printf FC "%s\n",qq(parmfile       $top   # Input PARM file);
printf FC "%s\n",qq(ambercoor      $rst   # Input coordinate file);
printf FC "%s\n",q(exclude        scaled1-4);
printf FC "%s\n",q(1-4scaling     0.833333  # =1/1.2, default is 1.0);
printf FC "%s\n",q();
printf FC "%s\n",q(# Follows the last step);
printf FC "%s\n",qq(velocities      $vel_pdb   ;# velocities from last run (pdb format));
printf FC "%s\n",qq(firsttimestep       $step              ;# last step of previous run);
printf FC "%s\n",q();
printf FC "%s\n",qq(outputname     $name # Prefix of output files);
printf FC "%s\n",q();
printf FC "%s\n",q(reinitvels     $temperature);
printf FC "%s\n",q();
printf FC "%s\n",q(run 500 ;# 1 ps (numsteps x timestep));
printf FC "%s\n",q();

#########################################################################
sub printHelpAndExit {
        printf STDERR "\n This script creates a NAMD configuration file to restart an AMBER MD run.\n\n";
	printf STDERR " Usage:\n";
	printf STDERR " $0 amber_topology_file amber_restart_file root_namd_files\n\n";
	exit;
};

#########################################################################
sub read_restart {
	my $f_in = shift @_;
	my $f_out = shift @_;
	open(F0,"$f_in") || die "Cannot open file $f_in: $!";
	my ($l1, $l2);
	chomp ($l1 = <F0>);
	chomp ($l2 = <F0>);
#line 2 contains number of atoms and step number:
	my ($natom, $step) = split(' ',$l2);
#number of lines is natom / 2 (6 coordinates per line)
	if ($natom % 2 != 0) { $natom++};
	my $nlines = $natom / 2;
#skip coordinates
	for my $i (1..$nlines){
		<F0>;
	}
#write velocites to file
	open(F1,">$f_out") || die "Cannot open file $f_out: $!";
	print F1 $l1,"\n";
	print F1 $l2,"\n";
	for my $i (1..$nlines){
		$l1 = <F0>;
		print F1 $l1;
	}
# get box dimensions from last line
	$l1 = <F0>;
	my ($b1, $b2, $b3, $ba1, $ba2, $ba3) = split(' ',$l1);
	close(F1);
	close(F0);
	return ($b1, $step);
}

#########################################################################
sub guess_pme {
	my $x = shift @_;
	my $int = int($x) + 4; # give some margin
	my $done = 0;
	while (! $done) {
		my @factors = prime_factors ($int);
##		foreach $i (@factors){print $i,"\n";};
		my @tmp = sort {$b<=>$a} @factors;
		if ($tmp[0] > 5) { $int++ }else{$done = 1};
	}
	return $int;
}

#########################################################################
sub prime_factors {
	my $number = shift @_;
	my $left = $number;
	my $test;
	my @list;
# loop through all possible factors
#
	foreach $test (2..$number){
 
# exit when no factoring left to do
#
		if ($left == 1) { last; }
 
# doest $test divide $left?
#
		if ($left % $test == 0) {
			$left /= $test;
			push @list, $test;
			redo ;
		}
	}
	return @list;
}

#########################################################################