Simulations

To run a simulation with SHARP Pack, only a single input file named param.in is required. This file contains all simulation details specified through predefined keywords.

Note

The file must be located in the same directory where the executable sharp.x is launched, unless a path is explicitly provided.

Keywords

Table Input Keywords lists all available input parameter keywords and their usage.

Table 7 Input Keywords
Keywords	Arguments	Descriptions
acckval	yes	If TRUE, writes accumulated result with k-value; useful for branching probability calculation(s) of Tully model
	\(\color{blue}{\text{no}}\)	No accumulated result(s)
approximation	\(\color{blue}{\text{CA}}\)	Centroid approximation RPSH (RPSH-CA)
	BA	Bead approximation RPSH (RPSH-BA)
dboc	\(\color{blue}{\text{BF}}\)	DBOC corrected force(s) at bead level
	CF	DBOC corrected force(s) at centroid level
decoherence	\(\color{blue}{\text{none}}\)	No decoherence correction
	damp	Expontial damping decoherence
dynamics	\(\color{blue}{\text{rpmd}}\)	RPMD dynamics with surface hopping
	\(\text{trpmd}\)	T-RPMD dynamics with surface hopping
iprint	\(n\)	Print data every n timesteps
iseed	\(n\)	Random number generator seed
nbead	\(n\)	Number of bead(s)
ncore	\(n\)	Number of core(s): (1-Serial), (n-Parallel)
nequil	\(n\)	equailibrate steps if pimd sampling
nmode	\(\color{blue}{\text{fft}}\)	FFT tranform for normal mode
	matrix	use normal mode tranformation matrix for normal mode
nparticle	\(n\)	Number of simulating particle(s)
nsample	\(n\)	sample run steps if pimd sampling
nstate	\(\color{blue}{n}\)	Number of electronic states; Default based on Model
nstep	\(n\)	Number of simulation steps
ntraj	\(n\)	Number of trajectories
model	\(s\)	Model name string from Table Models
pimd	\(\text{pile} \quad \tau\)	PIMD sampling of postion and velocity with PILE thermostat
pinit	\(f\)	Initial momentum (a.u.)
rsamp	fraction	Sample centroid bead from Gaussian distribution and other bead(s) position based on fraction of de-Broglie length
	\(\color{blue}{\text{gaussian}}\)	Sample initial bead(s) position from Gaussian distribution
	wigner	Sample initial bead(s) position from Wigner distribution
rmap	\(\color{blue}{\text{no}}\)	No mapping ring polymer centroid to dividing surface after sampling rsamp
	r0	Mapping ring polymer centroid to dividing surface after sampling rsamp
rmode	\(\color{blue}{\text{direct}}\)	Direct ring polymer sampling of rsamp
	norm	Normal mode ring polymer sampling of rsamp for gaussian
rundtail	yes	Print run-time details of the simulation
	\(\color{blue}{\text{no}}\)	No details of simulation
temperature	\(f\)	Temperature (K)
tstep	\(f\)	Simulation time-step (a.u.)
vsamp	fixed	Deterministic initial momentum
	\(\color{blue}{\text{gaussian}}\)	Sample initial momenta from Gaussian distribution
	wigner	Sample initial momenta from Wigner distribution
vrescale	\(\color{blue}{\text{CL}}\)	Centroid-approximation level of velocity rescaling of ring-polymer for energy conservation
	BL	Bead-approximation level of velocity rescaling of ring-polymer for energy conservation
vreverse	\(\color{blue}{\text{never}}\)	Never reversal (NR) velocity for frustrated hop(s)
	always	Always reversal (VR) of velocity for frustrated hop(s)
	delV1	Velocity reversal based on Truhlar’s \(\Delta V\) scheme
	delV2	Velocity reversal based on Subotnik’s \(\Delta V^2\) scheme
finish	***	End of the input file

Legend:: \(n\) – integer, \(f\) – real number, \(s\) – string, \(\color{blue}{\text{blue}}\) – denotes default value.

Table 8 Spin-boson model specific keywords
Keywords	Arguments	Descriptions
model	\(\text{spinboson} \quad \epsilon \quad \Delta \quad \text{enu}\)	Spin-boson model with parameters with scale factor `enu`
spectra	debye \(\quad E_r \quad \omega_c \quad T\)	Debye spectral density parameters and temperature
	ohmic \(\quad E_r \quad \omega_c \quad T\)	Ohmic spectral density parameters and temperature

Important

All input parameters are specified in atomic units (a.u.), except for temperature, which must be given in Kelvin.
For spin-boson parameters, values are expressed in energy unit of \(\text{cm}^{-1}\) including temperature, and scaled by the factor enu.
If enu = 1, the spin-boson parameters are interpreted directly as in atomic units.

Input File

A sample input file (param.in) for Spin-Boson model.

Listing 1 SAMPLE INPUT FILE (param.in)

 #SHARP Pack INPUT PARAMETERS
 #model(spinboson eps delta, enu)
 model  spinboson  0.0  1.0  104.25

 #bathspectrum (lambda wc temp)
 spectra  debye  1  1  2

 nbeads  8

 approx  CA

 nmode  matrix
 
 iseed  12345

 nParticle  100

 #number steps
 nsteps  41340

 #trajectories
 ntraj  1000

 ncore  20

 #timestep
 tstep  1.0

 #postion sampling 
 rsamp  gaussian
 rmode  direct
 rmap   no

 #velocity sampling
 vsamp  gaussian

 #velocity reversal scheme for frustrated hop
 vreverse  always

 #decoherence
 decoherence  none

 rundtail  yes

 finish
 #end of Input File

Output Files

After a successful simulation run, SHARP Pack generates the following output files:

Table 9 Output file(s)
File(s)	Descriptions
`param.out`	Lists all the model parameters used in the simulation.
`dcoupling.out`	Detailed trajectory data including positions, momenta, energies, NACs, etc.
`hoppinghist.out`	Records surface-hopping events between different potential energy surfaces.
`pop_adiabat1.out`	Adiabatic population computed using Method I.
`pop_adiabat2.out`	Adiabatic population computed using Method II.
`pop_diabat1.out`	Diabatic population computed using Method I.
`pop_diabat2.out`	Diabatic population computed using Method II.
`pop_diabat3.out`	Diabatic population computed using Method III.
`pop_branch.out`	Branching probablity of Tully model if `accKval` `yes` .
`sampling_pos.out`	Initial sampling distribution of positions.
`sampling_vel.out`	Initial sampling distribution of velocities.
`sampling_dis.out`	Initial sampling Boltzmann distribution function of speeds.
`pimd_sample.out`	PIMD sampling output (generated only if PIMD sampling is enabled).
`energy_surface.out`	Diabatic and adiabatic energy surfaces, including NAC vectors.
`energy_dboc.out`	Adiabatic and DBOC energy surfaces (if dboc is requested)
`psi.out`	Wavefunctions, adiabatic coefficients, and PES.
`bathfrequency.out`	Spectral density information (specific to the spin-boson model).

Note

The output files dcoupling.out, hoppinghist.out, and psi.out are produced only if rundtail is yes in the input file.

Population Calculation

SHARP Pack runs simulations in an adiabatic formalism; hence, the adiabatic state populations are the direct results. But it also provides multiple approaches for computing adiabatic and diabatic populations.

Adiabatic Population

There are two methods:

Method-I: Direct counting trajectory

The percentage of trajectories propagating on each adiabatic PES represents the population of that state.

\[P_{\alpha} = \langle \delta_{\alpha,\lambda} \rangle\]

where \(\lambda\) is the adiabatic active surface.

output file name: pop_adiabat1.out

Method-II: Density matrix method

Populations are computed from the complex-valued electronic coefficients.

\[P_{\alpha} = \langle \rho_{\alpha\alpha} \rangle\]

where \(\rho_{ij} = c_i\,c_j^*\) is the electronic density matrix.

output file name: pop_adiabat2.out

Diabatic Population

There are three methods, [21] using \(U_{ij}\) which are elements of the adiabatic-to-diabatic transformation matrix, and \(\lambda\) denotes the active surface.

Method-I: Projection Method

The diabatic population is determined by projecting the active adiabatic state onto the diabatic basis.

\[P_\alpha = \left\langle \sum_i |U_{\alpha i}\, \delta_{i,\lambda}|^2 \right\rangle,\]

output file name: pop_diabat1.out

Method-II: Electronic Wavefunction Method

Diabatic populations are calculated by transforming the electronic coefficients \(c_i\) in the adiabatic basis to the diabatic representation.

\[P_\alpha = \left\langle \sum_i | U_{\alpha i}\,c_i |^2 \right\rangle,\]

output file name: pop_diabat2.out

Method-III: Mixed quantum-classical density matrix

The diabatic populations are extracted from the electronic density matrix in the diabatic representation.

\[P_\alpha = \left\langle \sum_i |U_{\alpha i}|^2 \delta_{i,\lambda} + \sum_{i<j} 2\, \text{Re}\left[ U_{\alpha i} , \rho_{ij} , U_{\alpha j}^* \right] \right\rangle,\]

output file name: pop_diabat3.out

Running Simulation

To run SHARP Pack, a single input file param.in is required, which defines the system model parameters (see Input Keywords for details). Simulations are executed using the provided bash scripts located in the utility/ directory.

Available job submission scripts:

job-script-local.sh Run simulations on a local machine.
job-script-hpc.sh Submit jobs on an HPC cluster using the Slurm workload manager.
job-script-branching.sh Used only for computing branching probabilities on Tully model over multiple k-values on an HPC system.
job-script-avg.sh Submit job that computes average of output of the completed parallel jobs simulation.

Each script automatically generates the appropriate job submission file (serial or parallel, depending on the ncore value in param.in) and submits the job.

Usage

Run any of the scripts as follows:

# Run on a local machine
$ sh job-script-local.sh

# Run on an HPC cluster (Slurm)
$ sh job-script-hpc.sh

# Compute branching probability on HPC
$ sh job-script-branching.sh

Example script

An example of job_script_hpc.sh is shown below which is designed for running simulations on an HPC cluster using the Slurm workload manager. It automatically detects the number of cores (ncore) requested in the param.in file and prepares the submission accordingly:

If ncore 1; A serial job is created and submitted.
If ncore n (n > 1); The script creates n run* directories and distributes the number of trajectories across them for parallel simulation. After all simulations are complete, the results are averaged to produce the final output.

Listing 2 Bash script for job submission

#!/bin/bash
#============================================================
#     JOB SCRIPT TO RUN SHARP PACK v2 ON HPC SERIAL/PARALLEL
#     WHILE RUNNING PARALLEL, BASED ON n-CPU, IT CREATES
#     n-DIRECTORIES AND RUN JOBS PARALLELLY, AND FINALLY 
#     COMPUTES AVERAGE OF THE DISIRED POPULATATION (NEED TO
#     SPECIFY FOR INPUT, SEE EXAMPLE IN, WRITE_AVERAGE_INPUT)
#
#     authors    - D.K. Limbu & F.A. Shakib     
#
#     Method Development and Materials Simulation Laboratory
#     New Jersey Institute of Technology
#
#     USAGE:: bash job_script_hpc.sh
#           - creates and run job(s)
#
#============================================================

# Path to your executable
#============================================================
root2bin=~/Softwares/SHARP_pack_2/bin
exe=${root2bin}/sharp.x

#============================================================
maxcore=1
ncore=$(awk '/ncore/{print $2}' param.in)
maxcore=$((ncore > maxcore ? ncore : maxcore))
echo 'ncore(s):' $maxcore

# ==== CONFIGURATION ====
#============================================================
root_dir="run"            # Base directory name
ncore=$maxcore            # Number of parallel jobs/directories
param_file="param.in"     # Parameter file to copy to each directory

# Determine padding width based on ncore
if [ "$ncore" -lt 10 ]; then
   pad=1
elif [ "$ncore" -lt 100 ]; then
   pad=2
else
   pad=3
fi

# ==== FUNCTION TO CREATE AVERAGE INPUT ====
write_average_input() {
  #fname :: File name of specific population data for averaging
  #df    :: No. of data column(s)(except first column)?

  fname="pop_diabat3.out"   
  df=2

  echo "Calculating average results from $ncore parallel jobs."

  rm -f input

  echo ${root_dir} >> input
  echo ${fname} >> input
  echo ${ncore} >> input
  echo ${df} >> input
  echo ${pad} >> input

}

# ==== FUNCTION TO CREATE PARALLEL DIRECTORIES ====
paralleldir() {
  # Check if the first directory already exists
  first_dir=$(printf "${root_dir}%0${pad}g" 1)
  if [ -d "$first_dir" ]; then
    echo '**********************************'
    echo "$first_dir exists. Remove it first!"
    echo '**********************************'
    exit 1
  fi

  # Create directories and copy files
  for i in $(seq -f "%0${pad}g" 1 $ncore); do
    dir="${root_dir}${i}"
    mkdir "$dir"
    cp param.in "$dir/param.in"
    echo "$i" > "$dir/icpu.in"
  done

  echo "Created $ncore directories: ${root_dir}$(printf "%0${pad}d" 1) to ${root_dir}${ncore}"
}

# ==== FUNCTION TO WRITE SLURM ARRAY SCRIPT ====
write_array_job() {

  cat <<EOF > submit_array.sh
#!/bin/bash
#SBATCH --job-name=array_run
#SBATCH --array=1-${ncore}
#SBATCH --output=${root_dir}%0${pad}a/myoutput_%A_%a.out
#SBATCH --error=${root_dir}%0${pad}a/jobError_%A_%a.out
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=1
#SBATCH -p RM-shared
#SBATCH --mem-per-cpu=2000M
#SBATCH --mail-type=end,fail
#SBATCH --time=01:00:00

cd ${root_dir}\$(printf "%0${pad}d" \${SLURM_ARRAY_TASK_ID})

 $exe

EOF

  chmod +x submit_array.sh
}

# ==== FUNCTION TO WRITE SERIAL SLURM  SCRIPT ====
write_serial_job() {

  cat <<EOF > submit_serial.sh
#!/bin/bash
#SBATCH --job-name=serial_run
#SBATCH --output=myoutput%j.out
#SBATCH --error=jobError.out
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=1
#SBATCH -p RM-shared
#SBATCH --mem-per-cpu=2000M
#SBATCH --mail-type=end,fail
#SBATCH --time=10:00:00

 $exe

EOF

  chmod +x submit_serial.sh
}

# ==== FUNCTION TO WRITE ANALYSIS JOB SCRIPT ====
write_average_job() {
  cat <<EOF > average_job.sh
#!/bin/bash
#SBATCH --job-name=average
#SBATCH --output=myoutput%j.out
#SBATCH --error=jobError.out
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=1
#SBATCH -p RM-shared
#SBATCH --mem-per-cpu=2000M
#SBATCH --mail-type=end,fail
#SBATCH --time=00:10:00

  ${root2bin}/average.x

EOF

  chmod +x average_job.sh
}

#============================================================
## RUN SERIAL JOB ##
#============================================================
if [ $ncore -eq 1 ] ; then

  echo "Serial job running in a single core!!"

  write_serial_job

  sbatch submit_serial.sh

#============================================================
## RUN PARALLEL JOB(S) ON HPC ##
#============================================================
elif [ $ncore -gt 1 ] ; then
  #create parallel jobs

  echo "Parallel job running on " $ncore "core(s)."

  # ==== MAIN ====
  paralleldir
  write_array_job
  write_average_input
  write_average_job

  # Submit array job
  array_jobid=$(sbatch --parsable submit_array.sh)

  echo "Submitting post-processing job dependent on array job ${array_jobid}..."

  # Submit average job to run only after array job completes
  sbatch --dependency=afterok:$array_jobid average_job.sh

  echo "All jobs submitted successfully!"

fi

Note

To compute the average results from parallel runs, must to specify fname, and df inside the write_average_input() fucntion in job_script_hcp.sh or job_script_local.sh as:

# ==== FUNCTION TO CREATE AVERAGE INPUT ====
write_average_input() {
#fname :: File name of specific population data for averaging
#df    :: No. of data column(s)(except first column)

fname="pop_diabat3.out"
df=2
.
.

The final average output file is fname_ave.out. e.g. pop_diabat3.out ==> pop_diabat3_ave.out