#!usr/bin/perl
use strict;
require Exporter;
our @ISA = qw(Exporter);
our @EXPORT = qw(PSF, BasisVector, Minimize, MakeEquRef, Equilibrate, PULL );


#####################################################################################################################################
sub PSF
	{

my @imports = @_;

my $pdb = $imports[0];

chomp ($pdb);
open (OUT ,">make_psf.txt");

print OUT "
package require autopsf\n
mol new $pdb\.pdb\n
set A [atomselect top protein]
autopsf -mol 0 -protein -solvate -ionize\n
";
#~ autopsf -solvate -ionize -protein
system "vmd -dispdev text < make_psf.txt\n";
system "mv -v *_autopsf.pdb $pdb\_autopsf.pdb";
system "mv -v *_autopsf.psf $pdb\_autopsf.psf";
system "mv -v *_autopsf.log $pdb\_autopsf.log";
system "rm *autopsf_formatted*";
system "rm make_psf.txt\n";

my @break = split(/\//, $pdb);

	}
#####################################################################################################################################
sub BasisVector
	{
my @imports = @_;

my $pdb = $imports[0];

chomp ($pdb);


#ATOM   2514  OH2 TIP3W 328     -34.282  17.014   0.993  1.00  0.00      W1   O

system("grep TIP3W $pdb > water.txt");

open(WATER, "<water.txt");

my @x = ();
my @y = ();
my @z = ();

while(<WATER>)
	{
	my $in = $_;
	$in =~ s/TIP3W/  /;
	my @break = split(/ +/,$in);

	
	push @x,$break[4];
	push @y,$break[5];
	push @z,$break[6];
	
	
	}
my @x_sorted = sort {$a <=> $b} @x;
my @y_sorted = sort {$a <=> $b} @y;
my @z_sorted = sort {$a <=> $b} @z;
	
my $x_min = $x_sorted[0];
my $x_max = $x_sorted[-1];
my $y_min = $y_sorted[0];
my $y_max = $y_sorted[-1];
my $z_min = $z_sorted[0];
my $z_max = $z_sorted[-1];

my $x_length = $x_max - $x_min;
my $y_length = $y_max - $y_min;
my $z_length = $z_max - $z_min;
close(WATER);
return("$x_length", "$y_length", "$z_length");

	}
#####################################################################################################################################
sub Minimize
	{
my @imports = @_;

my $pdb = $imports[0];
my $cellBasisVector1 = $imports[1];
my $cellBasisVector2 = $imports[2];
my $cellBasisVector3 = $imports[3];
chomp ($pdb);
$pdb =~ s/_autopsf//g;

# Set up time for simulations.  MS for minimization, steps for MD.
	
my $MS = "10000";
my $steps = "0"; 
my $real_time = $steps/1000;
	
my $output = "$pdb\_minimize";
$output =~ s/\/pdb\/T0\d\d\d\//\/output\//g;
my $conf = $output;
$conf =~ s/\/output\//\/conf\//g;
open (MIN, ">$conf.conf");

# Generate the .conf file for minimization.

print MIN "


#############################################################
## JOB DESCRIPTION                                         ##
#############################################################

# This is what this job does
# N- C- Termini Constrant Velocity Pulling

#############################################################
## ADJUSTABLE PARAMETERS                                   ##
#############################################################

structure          $pdb\_autopsf.psf
coordinates        $pdb\_autopsf.pdb
outputName         $output

set temperature    50

# Continuing a job from the restart files
if {0} {
set inputname      $pdb
binCoordinates     \$inputname.restart.coor
binVelocities      \$inputname.restart.vel  ;# remove the \"temperature\" entry if you use this!
extendedSystem	   \$inputname.xsc
} 

firsttimestep      0


#############################################################
## SIMULATION PARAMETERS                                   ##
#############################################################

# Input parameter files.
paraTypeCharmm	    on
parameters          ./../models/par_all36_prot.prm 
parameters          ./../models/toppar_water_ions_namd.str
#parameters          ./../models/par_all36_carb.prm 
#parameters          ./../models/par_all36_lipid.prm 

# NOTE: Do not set the initial velocity temperature if you 
# have also specified a .vel restart file!
temperature         \$temperature
 

# Periodic Boundary conditions
# NOTE: Do not set the periodic cell basis if you have also 
# specified an .xsc restart file.
if {0} { 
cellBasisVector1    $cellBasisVector1  0    0
cellBasisVector2     0   $cellBasisVector2  0
cellBasisVector3     0    0   $cellBasisVector3
cellOrigin           0    0    0
}
wrapWater           on
wrapAll             on


# Force-Field Parameters
exclude             scaled1-4
1-4scaling          1.0
cutoff              12.0
switching           on
switchdist          10.0
pairlistdist        14.0


# Integrator Parameters
timestep            1.0  ;# 1fs/step
rigidBonds          water  ;# TIP3P
nonbondedFreq       2
fullElectFrequency  4  
stepspercycle       8

# Constant Temperature Control
langevin            on    ;# do langevin dynamics
langevinDamping     1     ;# damping coefficient (gamma) of 5/ps
langevinTemp        \$temperature
langevinHydrogen    no    ;# don't couple langevin bath to hydrogens

# Constant Pressure Control (variable volume)
if {1} {
useGroupPressure      no ;# needed for 2fs steps
useFlexibleCell       no  ;# no for water box, yes for membrane
useConstantArea       no  ;# no for water box, yes for membrane

langevinPiston        on
langevinPistonTarget  1.01325 ;#  in bar -> 1 atm
langevinPistonPeriod  100.0
langevinPistonDecay   50.0
langevinPistonTemp    \$temperature
}

#Output frequences for restart and output files.
restartfreq         1000     ;# 1000steps = every 1ps
dcdfreq             1000
xstFreq             1000
outputEnergies      1000
outputPressure      1000

#############################################################
## EXECUTION SCRIPT                                        ##
#############################################################

# Minimization
if {1} {
minimize            $MS
reinitvels          \$temperature
}

run $steps ;# $real_time ps

";


print "minimization run:\n"; 
print "./namd2 $conf.conf > $conf.log\n";

#~ system("./namd2 $conf.conf > $conf.log\n");

#~ open(UPDATE, ">run_min.txt");
#~ print UPDATE "mol new $pdb\_autopsf.psf\n\nmol addfile $pdb\_autopsf.pdb\n\nmol addfile $pdb\_minimize.dcd waitfor all\n\nset A [atomselect top all]\n\n\$A writepdb $pdb\_minimize_autopsf.pdb\n\n";
#~ close(UPDATE);

#~ system("vmd -dispdev text < run_min.txt");

	}
#####################################################################################################################################
sub MakeEquRef
{
	my @inputs = @_;
	my $pdb = $inputs[0];
	chomp($pdb);
	$pdb =~ s/_autopsf//g;
	open(IN, "< $pdb\_autopsf.pdb");
	open(OUT, "> $pdb\_equ.ref");
	while(<IN>)
	{
		my $line = $_;
		if ($line =~ m/ CA /g)
		{
			$line =~ s/1.00  0.00/0.00  1.00/g;
			$line =~ s/0.00  0.00/0.00  1.00/g;
		}
		else
		{
			$line =~ s/1.00  1.00/0.00  0.00/g;
			$line =~ s/1.00  0.00/0.00  0.00/g;
			$line =~ s/0.00  1.00/0.00  0.00/g;
		}
		print OUT "$line";
	}
}
#####################################################################################################################################
sub Equilibrate
	{	
my @inputs = @_;

my $MS = $inputs[0];
my $steps = $inputs[1];
my $real_time = $inputs[2];
my $pdb = $inputs[3];
my $cellBasisVector1 = $inputs[4];
my $cellBasisVector2 = $inputs[5];
my $cellBasisVector3 = $inputs[6];
chomp($pdb);
$pdb =~ s/_autopsf//g;

my $MS = 0;
my $steps = 500000; 
$real_time = $steps/1000;
my $output = "$pdb\_Equilibrate";
$output =~ s/\/pdb\/T0\d\d\d\//\/output\//g;
my $input = $output;
$input =~ s/Equilibrate/minimize/g;
my $conf = $output;
$conf =~ s/\/output\//\/conf\//g;

open(EQU, ">$conf.conf");

print EQU "
#############################################################
## JOB DESCRIPTION                                         ##
#############################################################

# This is what this job does
# N- C- Termini Constrant Velocity Pulling

#############################################################
## ADJUSTABLE PARAMETERS                                   ##
#############################################################

structure          $pdb\_autopsf.psf
coordinates        $pdb\_autopsf.pdb
outputName         $output

set temperature    10

# Continuing a job from the restart files
if {1} {
set inputname      $input
binCoordinates     \$inputname.restart.coor
binVelocities      \$inputname.restart.vel  ;# remove the \"temperature\" entry if you use this!
extendedSystem	   \$inputname.restart.xsc
} 

firsttimestep      0


#############################################################
## SIMULATION PARAMETERS                                   ##
#############################################################

# Input
paraTypeCharmm	    on
parameters          ./../models/par_all36_prot.prm 
parameters          ./../models/toppar_water_ions_namd.str
#parameters          ./../models/par_all36_carb.prm 
#parameters          ./../models/par_all36_lipid.prm 

# NOTE: Do not set the intial velocity temperature if you 
# have also specified a .vel restart file.
#~ temperature         \$temperature
 

# Periodic Boundary conditions
# NOTE: Do not set the periodic cell basis if you have also 
# specified an .xsc restart file!
if {0} { 
cellBasisVector1    $cellBasisVector1  0    0
cellBasisVector2     0   $cellBasisVector1  0
cellBasisVector3     0    0   $cellBasisVector1
cellOrigin           0    0    0
}
wrapWater           on
wrapAll             on


# Force-Field Parameters
gbis on
exclude             scaled1-4
1-4scaling          1.0
cutoff              12.0
switching           on
switchdist          10.0
pairlistdist        14.0


# Integrator Parameters
timestep            1.0  ;# 1fs/step
rigidBonds          water  ;# TIP3P
nonbondedFreq       2
fullElectFrequency  4  
stepspercycle       8

# Temperature Parateters
reassignFreq 1000
reassignIncr 10 
reassignTemp 10
reassignHold 300

#PME (for full-system periodic electrostatics)
if {0} {
PME                 yes
PMEGridSpacing      1.0

#manual grid definition
#PMEGridSizeX        32
#PMEGridSizeY        32
#PMEGridSizeZ        64
}

# Constant Pressure Control (variable volume)
if {1} {
useGroupPressure      no  ;# needed for 2fs steps
useFlexibleCell       no  ;# no for water box, yes for membrane
useConstantArea       no  ;# no for water box, yes for membrane

langevinPiston        on
langevinPistonTarget  1.01325 ;#  in bar -> 1 atm
langevinPistonPeriod  100.0
langevinPistonDecay   50.0
langevinPistonTemp    300
}

restartfreq         1000     ;# 1000steps = every 1ps
dcdfreq             1000
xstFreq             1000
outputEnergies      1000
outputPressure      1000

constraints on
conskfile $pdb\_equ.ref
consref $pdb\_equ.ref
conskcol B

#############################################################
## EXECUTION SCRIPT                                        ##
#############################################################

#~ minimize            1000
run $steps ;# $real_time ps

";
print "./namd2 $conf\_equ.conf > $conf\_equ.log\n";
#system("./namd2 $conf\_equ.conf > $conf\_equ.log\n");

#set up for step2 restart simulations in case of cell crash.
my $output = "$pdb\_Equilibrate_step2";
$output =~ s/\/pdb\/T0\d\d\d\//\/output\//g;
$input =~ s/minimize/Equilibrate/g;

close(EQU);

open(EQU, "> $conf\_step2.conf");
print EQU "
#############################################################
## JOB DESCRIPTION                                         ##
#############################################################

# This is what this job does
# N- C- Termini Constrant Velocity Pulling

#############################################################
## ADJUSTABLE PARAMETERS                                   ##
#############################################################

structure          $pdb\_autopsf.psf
coordinates        $pdb\_autopsf.pdb
outputName         $output

#set temperature    10

set inputname      $pdb

# Continuing a job from the restart files
if {1} {
binCoordinates     $input.restart.coor
binVelocities      $input.restart.vel  ;# remove the \"temperature\" entry if you use this!
extendedSystem	   $input.restart.xsc
}
# procedure to get the first time step for the new simulation
# from the old simulation whether it be a completed sim or one
# to be restarted
if {0} {
proc get_first_ts { xscfile } {
  set fd [open \$xscfile r]
  gets \$fd
  gets \$fd
  gets \$fd line
  set ts [lindex \$line 0]
  close \$fd
  return \$ts
}
if [file exists \$inputname.xsc] {
set firsttime [get_first_ts \$inputname.xsc]
firsttimestep \$firsttime
} else {
set firsttime [get_first_ts \$inputname.restart.xsc]
firsttimestep \$firsttime
}

} else {
firsttimestep      0
}


#############################################################
## SIMULATION PARAMETERS                                   ##
#############################################################

# Input
paraTypeCharmm	    on
parameters          ./../models/par_all36_prot.prm 
parameters          ./../models/toppar_water_ions_namd.str
#parameters          ./../models/par_all36_carb.prm 
#parameters          ./../models/par_all36_lipid.prm 

# NOTE: Do not set the initial velocity temperature if you 
# have also specified a .vel restart file!
#temperature         \$temperature
 

# Periodic Boundary conditions
# NOTE: Do not set the periodic cell basis if you have also 
# specified an .xsc restart file!
if {0} { 
cellBasisVector1    $cellBasisVector1  0    0
cellBasisVector2     0   $cellBasisVector1  0
cellBasisVector3     0    0   $cellBasisVector1
cellOrigin           0    0    0
}
wrapWater           on
wrapAll             on


# Force-Field Parameters
gbis on
exclude             scaled1-4
1-4scaling          1.0
cutoff              12.0
switching           on
switchdist          10.0
pairlistdist        14.0


# Integrator Parameters
timestep            1.0  ;# 1fs/step
rigidBonds          water  ;# TIP3P
nonbondedFreq       2
fullElectFrequency  4  
stepspercycle       8

# Temperature Parateters
reassignFreq 1000
reassignIncr 10 
reassignTemp 10
reassignHold 300

#PME (for full-system periodic electrostatics)
if {0} {
PME                 yes
PMEGridSpacing      1.0

#manual grid definition
#PMEGridSizeX        32
#PMEGridSizeY        32
#PMEGridSizeZ        64
}

# Constant Pressure Control (variable volume)
if {0} {
useGroupPressure      no ;# needed for 2fs steps
useFlexibleCell       no  ;# no for water box, yes for membrane
useConstantArea       no  ;# no for water box, yes for membrane

langevinPiston        on
langevinPistonTarget  1.01325 ;#  in bar -> 1 atm
langevinPistonPeriod  100.0
langevinPistonDecay   50.0
langevinPistonTemp    300
}

restartfreq         1000     ;# 1000steps = every 1ps
dcdfreq             1000
xstFreq             1000
outputEnergies      1000
outputPressure      1000

constraints on
conskfile $pdb\_equ.ref
consref $pdb\_equ.ref
conskcol B

#############################################################
## EXECUTION SCRIPT                                        ##
#############################################################

run $steps ;# $real_time ps

";
print "RESTART RUN FROM CRASH\n";
print "./namd2 $pdb\_equ.conf > $pdb\_equ_step2.log\n";
#system("./namd2 $pdb\_equ.conf > $pdb\_equ_step2.log\n");

	}
#####################################################################################################################################
sub PULL
	{

my @inputs = @_;

my $MS = $inputs[0];
my $steps = $inputs[1];
my $real_time = $inputs[2];
my $pdb = $inputs[3];
chomp($pdb);
my $cellBasisVector1 = $inputs[4];
my $cellBasisVector2 = $inputs[5];
my $cellBasisVector3 = $inputs[6];

open (OUT ,">make_pdb.txt");
open (UPDATE, "> min.txt");
open(PDB, "<$pdb\_autopsf.pdb");
my $CA_flag = 1;
open(REF, ">$pdb.ref");
my $last_CA;

while(<PDB>)
	{
	my $in = $_;
	my @break = split(/ +/, $in);
	my $atom_number = $break[1];
	my $atom_type = $break[2];
	if($atom_type eq "CA")
		{
		$last_CA = $atom_number;
		}
	}
	
close(PDB);
open(PDB, "<$pdb\_autopsf.pdb");

my $line;
open(REF, ">$pdb.ref");
while(<PDB>)
	{
	my $in = $_;
	chomp($in);
	my $last_res_number;
	my @break = split(/1.00  0.00/, $in);
	my $First = $break[0];
	my $Second = $break[1];
	my $occupancy;
	my $B;
	$B = "0.00";
	$occupancy = "0.00";
	if($First =~ m/TIP3W/g || $Second =~ m/ION/g)
		{
		print REF "$First$occupancy  $B$Second\n"
		}
	else
		{

		if($First =~ m/ATOM/g)
			{
	
			$occupancy = "0.00";
			
			if($First =~ m/CA/g && $CA_flag == 1)
				{
				$B = "1.00";
				$CA_flag = 0;
				}
			else
				{
				$B = "0.00";
				}
			if($First =~ s/$last_CA  CA/$last_CA  CA/g)
				{
				$occupancy = "1.00";
				}
			$line = "$First$occupancy  $B$Second";
			$line =~ s/0\.00  0\.00ATOM/ATOM/g;
			print REF "$line\n";
			}
		}
	}
close(REF);

print "$pdb.ref\n";
open(REF, "<$pdb.ref");
my $x1;
my $y1;
my $z1;
my $x2;
my $y2;
my $z2;

while(<REF>)
	{
	my $in = $_;
	my $last_res_number;
	my @break = split(/ +/, $in);
	my $atom  = $break[0];
	my $atom_number = $break[1];
	my $atom_type = $break[2];
	my $res_type = $break[3];
	my $chain = $break[4];
	my $res_number = $break[5];
	my $x = $break[6];
	my $y = $break[7];
	my $z = $break[8];
	my $occupancy = $break[9];
	my $B = $break[10];
	my $other = $break[11];
	
	# get one end
	if($B == 1.00)	
		{
		$x1 = $x;
		$y1 = $y;
		$z1 = $z;
		}
	# get other end
	if($occupancy == 1.00)
		{
		$x2 = $x;
		$y2 = $y;
		$z2 = $z;
		}
	}

# calculate pulling vector	
my $dx = $x2 - $x1; my $dy = $y2 - $y1; my $dz = $z2 - $z1;
my $dx_norm; my $dy_norm; my $dz_norm;
my $s = sqrt(($dx * $dx) + ($dy * $dy) + ($dz * $dz));

if($s > 0)
	{	
	$dx_norm = $dx / $s;
	$dy_norm = $dy / $s;
	$dz_norm = $dz / $s;
	}
else
	{	
	$dx_norm = 1;
	$dy_norm = 1;
	$dz_norm = 1;
	}

#$SMDK is the spring constant
#$SMDVel is the velocity of the SMD atoms
#$SMDX is the x component
#$SMDY is the y component
#$SMDZ is the z component

my $SMDK = 7;
my $SMDVel = 0.001;
my $SMDX = $dx_norm;
my $SMDY = $dy_norm;
my $SMDZ = $dz_norm;

my $output = "$pdb\_pulled_c_v";
$output =~ s/\/pdb\/T0\d\d\d\//\/output\//g;
my $input = $output;
$input =~ s/pulled_c_v/Equilibrate_step2/g;
my $conf = $output;
$conf =~ s/\/output\//\/conf\//g;

open(NAMD, ">$conf.conf");

my @break_fix_for_pdb_name = split(/\//, $pdb);
my $name = $break_fix_for_pdb_name[2];
print NAMD "

#############################################################
## JOB DESCRIPTION                                         ##
#############################################################

# This is what this job does
# N- C- Termini Constrant Velocity Pulling

#############################################################
## ADJUSTABLE PARAMETERS                                   ##
#############################################################

structure          $pdb\_autopsf.psf
coordinates        $pdb\_autopsf.pdb
outputName         $output

#set temperature    300

# Continuing a job from the restart files
if {1} {
set inputname      $input
binCoordinates     \$inputname.restart.coor
binVelocities      \$inputname.restart.vel  ;# remove the \"temperature\" entry if you use this!
extendedSystem	   \$inputname.restart.xsc
} 

firsttimestep      0


#############################################################
## SIMULATION PARAMETERS                                   ##
#############################################################

# Input
paraTypeCharmm	    on
parameters          ./../models/par_all36_prot.prm 
parameters          ./../models/toppar_water_ions_namd.str
#parameters          ./../models/par_all36_carb.prm 
#parameters          ./../models/par_all36_lipid.prm 

# NOTE: Do not set the initial velocity temperature if you 
# have also specified a .vel restart file!
#temperature         \$temperature
 

# Periodic Boundary conditions
# NOTE: Do not set the periodic cell basis if you have also 
# specified an .xsc restart file!
if {0} { 
cellBasisVector1    $cellBasisVector1  0    0
cellBasisVector2     0   $cellBasisVector2  0
cellBasisVector3     0    0   $cellBasisVector3
cellOrigin           0    0    0
}
wrapWater           on
wrapAll             on


# Force-Field Parameters
gbis on
exclude             scaled1-4
1-4scaling          1.0
cutoff              12.0
switching           on
switchdist          10.0
pairlistdist        14.0


# Integrator Parameters
timestep            1.0  ;# 1fs/step
rigidBonds          water  ;# TIP3P
nonbondedFreq       2
fullElectFrequency  4  
stepspercycle       8


#PME (for full-system periodic electrostatics)
if {0} {
PME                 yes
PMEGridSpacing      1.0

#manual grid definition
#PMEGridSizeX        32
#PMEGridSizeY        32
#PMEGridSizeZ        64
}


# Constant Temperature Control
langevin            on    ;# do langevin dynamics
langevinDamping     1     ;# damping coefficient (gamma) of 5/ps
langevinTemp        300
langevinHydrogen    no    ;# don't couple langevin bath to hydrogens


# Constant Pressure Control (variable volume)
if {1} {
useGroupPressure      no ;# needed for 2fs steps
useFlexibleCell       no  ;# no for water box, yes for membrane
useConstantArea       no  ;# no for water box, yes for membrane

langevinPiston        on
langevinPistonTarget  1.01325 ;#  in bar -> 1 atm
langevinPistonPeriod  100.0
langevinPistonDecay   50.0
langevinPistonTemp    300
}


restartfreq         1000     ;# 1000steps = every 1ps
dcdfreq             1000
xstFreq             1000
outputEnergies      1000
outputPressure      1000


# Fixed Atoms Constraint (set PDB beta-column to 1)
if {1} {
fixedAtoms          on
fixedAtomsFile      $pdb.ref
fixedAtomsCol       B
}





#############################################################
## EXTRA PARAMETERS                                        ##
#############################################################

# Put here any custom parameters that are specific to 
# this job (e.g., SMD, TclForces, etc...)
SMD on
SMDFile $pdb.ref
SMDk $SMDK
SMDVel $SMDVel
SMDDir $SMDX $SMDY $SMDZ
SMDOutputFreq 8


# to help produce something for the next run if this one crashes
###restartfreq <not sure best value for this>

#############################################################
## EXECUTION SCRIPT                                        ##
#############################################################

# Minimization
if {0} {
minimize            $MS
#reinitvels          \$temperature
}

run $steps ;# $real_time ps

";


#print "SMD run:\n";
print "./namd2 $conf\_pulled_c_v.conf > $conf\_pulled_c_v.log\n";
#system("./namd2 $pdb\_pulled_c_v.conf > $pdb\_pulled_c_v.log\n");

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


#my $output = "$pdb\_pulled_c_v_step2";

$output =~ s/_pulled_c_v/_pulled_c_v_step2/;
$input =~ s/Equilibrate_step2/pulled_c_v/g;
$conf =~ s/_pulled_c_v/_pulled_c_v_step2/;

close(NAMD);

open(NAMD, ">$conf.conf");

print NAMD "

#############################################################
## JOB DESCRIPTION                                         ##
#############################################################

# This is what this job does
# N- C- Termini Constrant Velocity Pulling

#############################################################
## ADJUSTABLE PARAMETERS                                   ##
#############################################################

structure          $pdb\_autopsf.psf
coordinates        $pdb\_autopsf.pdb
outputName         $output

#set temperature    300

set inputname      $pdb

# Continuing a job from the restart files
if {1} {
binCoordinates     $input.restart.coor
binVelocities      $input.restart.vel  ;# remove the \"temperature\" entry if you use this!
extendedSystem	   $input.restart.xsc
}
# procedure to get the first time step for the new simulation
# from the old simulation whether it be a completed sim or one
# to be restarted
if {0} {
proc get_first_ts { xscfile } {
  set fd [open \$xscfile r]
  gets \$fd
  gets \$fd
  gets \$fd line
  set ts [lindex \$line 0]
  close \$fd
  return \$ts
}
if [file exists \$inputname.xsc] {
set firsttime [get_first_ts \$inputname.xsc]
firsttimestep \$firsttime
} else {
set firsttime [get_first_ts \$inputname.restart.xsc]
firsttimestep \$firsttime
}

} else {
firsttimestep      0
}


#############################################################
## SIMULATION PARAMETERS                                   ##
#############################################################

# Input
paraTypeCharmm	    on
parameters          ./../models/par_all36_prot.prm 
parameters          ./../models/toppar_water_ions_namd.str
#parameters          ./../models/par_all36_carb.prm 
#parameters          ./../models/par_all36_lipid.prm 

# NOTE: Do not set the initial velocity temperature if you 
# have also specified a .vel restart file!
#temperature         \$temperature
 

# Periodic Boundary conditions
# NOTE: Do not set the periodic cell basis if you have also 
# specified an .xsc restart file!
if {0} { 
cellBasisVector1    $cellBasisVector1  0    0
cellBasisVector2     0   $cellBasisVector2  0
cellBasisVector3     0    0   $cellBasisVector3
cellOrigin           0    0    0
}
wrapWater           on
wrapAll             on


# Force-Field Parameters
gbis on
exclude             scaled1-4
1-4scaling          1.0
cutoff              12.0
switching           on
switchdist          10.0
pairlistdist        14.0


# Integrator Parameters
timestep            1.0  ;# 1fs/step
rigidBonds          water  ;# TIP3P
nonbondedFreq       2
fullElectFrequency  4  
stepspercycle       8


#PME (for full-system periodic electrostatics)
if {0} {
PME                 yes
PMEGridSpacing      1.0

#manual grid definition
#PMEGridSizeX        32
#PMEGridSizeY        32
#PMEGridSizeZ        64
}




# Constant Pressure Control (variable volume)
if {0} {
useGroupPressure      no ;# needed for 2fs steps
useFlexibleCell       no  ;# no for water box, yes for membrane
useConstantArea       no  ;# no for water box, yes for membrane

langevinPiston        on
langevinPistonTarget  1.01325 ;#  in bar -> 1 atm
langevinPistonPeriod  100.0
langevinPistonDecay   50.0
langevinPistonTemp    300
}


restartfreq         1000     ;# 1000steps = every 1ps
dcdfreq             1000
xstFreq             1000
outputEnergies      100
outputPressure      1000


# Fixed Atoms Constraint (set PDB beta-column to 1)
if {1} {
fixedAtoms          on
fixedAtomsFile      $pdb.ref
fixedAtomsCol       B
}


# IMD Settings (can view sim in VMD)
if {0} {
IMDon           on
IMDport         3000    ;# port number (enter it in VMD)
IMDfreq         1       ;# send every 1 frame
IMDwait         no      ;# wait for VMD to connect before running?
}


#############################################################
## EXTRA PARAMETERS                                        ##
#############################################################

# Put here any custom parameters that are specific to 
# this job (e.g., SMD, TclForces, etc...)
SMD on
SMDFile $pdb.ref
SMDk $SMDK
SMDVel $SMDVel
SMDDir $SMDX $SMDY $SMDZ
SMDOutputFreq 8


#############################################################
## EXECUTION SCRIPT                                        ##
#############################################################



run $steps ;# $real_time ps

";

print "RESTART FROM CRASH SMD run:\n";
print "./namd2 $conf.conf > $conf.log\n";
#system("./namd2 $pdb\_pulled_c_v_step2.conf > $pdb\_pulled_c_v_step2.log\n");
	
	}
#####################################################################################################################################

1;