####################################################
# at T (K) and initial configuration with known density corresponding pressure 
#################### SETTING UP ####################
units real
dimension 3
newton on
boundary p p p
atom_style full

#################### VARIABLES ####################
variable        ext_temp equal ${temp}
variable MolMass equal 44.01    # Average molar mass (kg/kmol)
variable MolAtom equal 3        # Number of atoms in a molecule

variable tstep equal 1.0        # 1fs 
variable Ninptpre equal 10000   # Pre-initialize the NPT ensemble (very small timesteps)
variable Ninpt equal 50000      # Initialize the NPT ensemble
variable Npnpt equal 2000000    # Production in the NPT ensemble (volume)
variable Ninvtpre equal 10000   # Pre-initialize the NVT ensemble (very small timesteps)
variable Ninvt equal 50000      # Initialize the NVT ensemble
variable Npnvt equal 10000000    # Production in the NVT ensemble (energy)
variable Nrun equal 1000000     # production in the NVE ensemble
variable Nf equal 10000         # Nfreq (fix ave/time and thermo)
variable Ne equal  10           # Nevery (fix ave/time)
variable Nr equal ${Nf}/${Ne}   # Nrepeat (fix ave/time)
variable NBR equal ${Npnpt}/10  # Block averaging for density (NPT)
variable Nd equal ${Npnvt}/10   # Frequency of outputting positions of atoms in the NVT ensemble

variable nb equal 10            # Number of blocks (fix ordern)
variable nbe equal 20           # Number of block elements (fix ordern)
variable Nvisc equal 5          # Nevery for viscosity (fix ordern: sample data at (Nvisc*2))
variable Ncond equal 5          # Nevery for Tconductivity (fix ordern: sample data at (Ncond*2))
variable Ndiff equal 1000       # Nevery for diffusivity (fix ordern: sample data at (Ndiff))
variable Nwrit equal 100000     # Nwrite for transport properties (fix ordern: write every (Nwrit))

#################### ATOM DEFINITION and FORCEFIELD ####################
read_data       ${infile}
replicate	1 1 1				     # cell replication (none = 1 1 1)
include    forcefieldTraPPE.data                     # read the force field
delete_bonds     all bond 1 
delete_bonds     all angle 1 
#################### INITIALIZATION ####################
# neighbor lists
#neighbor 2.0 bin
#neigh_modify every 1 delay 0 check no 

# initializing velocities
velocity all create ${ext_temp} 432567 dist uniform

reset_timestep 0

# rate of writing thermal properties to the log file
thermo ${Nf}

#################### 3) Initializing the NPT ensemble ####################
# applying the shake or rigid algorithm for rigid molecules and the NPT solver
#fix constrain all shake 1.0e-6 1000 0 b 1 

fix integrate all rigid/npt molecule temp ${ext_temp} ${ext_temp} 10.0 iso 94.0 94.0 100.0

# Initializing the whole system with very small timesteps in the NPT ensemble
timestep 0.001
run ${Ninptpre} 
timestep 0.01
run ${Ninptpre} 
timestep 0.1
run ${Ninptpre} 
timestep 0.2
run ${Ninptpre} 
timestep 0.5
run ${Ninptpre} 
reset_timestep 0

# continuing the initialization with the final value of timestep
unfix integrate

fix integrate all rigid/npt molecule temp ${ext_temp} ${ext_temp} 100.0 iso 94.0 94.0 1000.0

timestep ${tstep}
run ${Ninpt}
reset_timestep 0

#################### 2) Obtaining average volume in NPT ######################
# Getting the average volume of the system
variable Volume equal vol
fix VoluAve all ave/time 1 ${Npnpt} ${Npnpt} v_Volume file volume.dat

# Getting the average density of the system (block averaging)
variable nAvog equal 6.0221415e23                     # Avogadro's number
variable A3_in_m3 equal 1e-30                         # Angstrom^3 in m^3
variable nMolecule equal atoms/${MolAtom}             # Total number of molecules
variable mMolecule equal (${MolMass}*1e-3/${nAvog})   # mass of a molecule (kg/molecule)
variable Dens equal v_nMolecule*${mMolecule}/(vol*${A3_in_m3})
fix DensAve all ave/time 1 ${NBR} ${NBR} v_Dens file density.dat

thermo          100                                    # Specify the interval between screen output of thermodynamic averages
thermo_style    custom step temp density epair press   # Format for screen output of thermodynamics

dump 	     1 all xyz 200000 dump_npt.xyz
dump_modify  1  element C1 O2 

run ${Npnpt}

quit
#################### 3) Initializing the NVT ensemble ####################
# scaling the size of the system to the average volume
variable sidesize equal (f_VoluAve^(1.0/3.0))    # get the volume
variable xlow equal xlo
variable ylow equal ylo
variable zlow equal zlo
variable xhig equal (xlo+${sidesize})
variable yhig equal (ylo+${sidesize})
variable zhig equal (zlo+${sidesize})
change_box all x final ${xlow} ${xhig} y final ${ylow} ${yhig} z final ${zlow} ${zhig}
unfix DensAve
unfix VoluAve

# changing the ensemble to NVT

undump 1
reset_timestep 0

unfix integrate
fix integrate all rigid/nvt molecule temp ${ext_temp} ${ext_temp} 100.0

# Initializing the whole system with very small timesteps in the NVT ensemble
timestep 0.001
run ${Ninvtpre} 
timestep 0.01
run ${Ninvtpre} 
timestep 0.1
run ${Ninvtpre} 
timestep 0.2
run ${Ninvtpre} 
timestep 0.5
run ${Ninvtpre} 
reset_timestep 0

# continuing the initialization with the final value of timestep
timestep ${tstep}
run ${Ninvt}
reset_timestep 0
#################### 4) Obtaining average total energy in NVT and the PVT data #####################
variable T1 equal temp
variable TE1 equal etotal
variable KE1 equal ke
variable PE1 equal pe
fix Tave1 all ave/time ${Npnvt} 1 ${Npnvt} v_T1
fix TEave1 all ave/time 1 ${Npnvt} ${Npnvt} v_TE1
fix KEave1 all ave/time ${Npnvt} 1 ${Npnvt} v_KE1
fix PEave1 all ave/time ${Npnvt} 1 ${Npnvt} v_PE1

###Output Style

variable        Nevery equal 20                      #Specify the distance (in timesteps) between samples for computing ensemble aver$
                                                     #  (Interval between samples in a block)
variable        Nrepeat equal 50                     #Specify the number of samples per output  of thermodynamic averages
                                                     #  (Samples per block output)
variable        Nfreq equal ${Nevery}*${Nrepeat}     #Specify the dump interval (in timesteps)
                                                     #  (Some people call these blocks)
variable        PotEpair equal epair          #Call the intermolecular potential energy $PotentialEnergy
variable        TotalEnergy equal etotal
variable        Pressure equal press                 #Call the pressure $Pressure
variable        Temperature equal temp               #Call the instantaneous temperature $Temperature
variable        Density equal density                #Call the instantaneous density $Density (will not change for NVT)

fix             2 all ave/time ${Nevery} ${Nrepeat} ${Nfreq} v_Temperature v_Density v_PotEpair v_TotalEnergy v_Pressure file ave.dens18.0molL.dat

thermo          100                                    # Specify the interval between screen output of thermodynamic averages
thermo_style    custom step temp density epair press   # Format for screen output of thermodynamics
                                                       # (MD step, temperature, internal energy, pressure) 
dump 	     2 all xyz 500000 dump_nvt.xyz
dump_modify  2  element C1 O2  
run ${Npnvt}