Simulation of THz-driven processes using LAMMPS

©️ Copyright 2025 @ Yao Wu
Author： Yao Wu 📨
Date：2025-10-20
Lisence：This document is licensed under Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0) license.

1 Introduction

This tutorial is to give a guideline to calculate THz-driven processes by using Lammps.

Terahertz (THz, 0.1~10 THz ) photons carry only a few millielectronvolts of energy (4.1 meV at 1 THz)—far below the band gaps, lattice-vibration phonon energies, or chemical-bond strengths of most materials. Consequently, THz light does not directly drive inter-band electronic transitions, yet it can resonantly couple to low-energy degrees of freedom such as phonons, free carriers, magnons, polarons, and superconducting gaps, thereby “gently” steering collective behavior far from equilibrium. In recent years, the peak power of THz sources (synchrotron radiation, ultrafast lasers, free-electron lasers, and intense-field THz generators) has surpassed GW cm⁻², giving rise to a wealth of novel condensed-matter phenomena.

The lammps simulation of terahertz does not directly introduce a light field. Instead, it adds an electric field in the shape of a THz Gaussian packet on the basis of the ordinary simulation. Since this band rarely involves electron transitions, only the ions can be applied with an electric field.

Below, we will take the simple BaTiO₃ (BTO) system as a test case and investigate its THz-driven processes.

2 Data preparation

Born effective charge calculation

Before applying the THz field, we need to compute the Born effective charges (Z* ) of the given system from DFT. For every ion, Z* is obtained by summing the appropriate tensor elements, weighted by the eigen-vector of the IR-active mode; this guarantees that the anions and cations carry equal and opposite charges, so that the unit cell remains neutral overall.

The basic input parameters in VASP are as follows, and those in other first-principles calculation software should be similar.

SYSTEM=BTO
ISTART=0
EDIFF = 1E-8
LWAVE=F
LCHARG=F
ENCUT = 500
LEPSILON = .TRUE.
LREAL = .FALSE.
NSW = 0
PREC = A

For simple systems such as BTO, the Born effective charge can be extracted through the following script:

#!/bin/bash

# Define input and output files
input_file="OUTCAR"
output_file="born.dat"

# Extract the BORN EFFECTIVE CHARGES block with grep and clean it with sed
grep -A 10000 "BORN EFFECTIVE CHARGES" "$input_file" | \
sed '/^$/d; /^[[:alpha:]]/d; /ion/d' > "$output_file"

# Re-format born.dat: keep only the first four columns for atoms 1, 2, 3
awk '{if ($1 == 1) print $1, $2, $3, $4}' "$output_file" > "formatted_born.dat"
echo -e "\n" >> "formatted_born.dat"
awk '{if ($1 == 2) print $1, $2, $3, $4}' "$output_file" >> "formatted_born.dat"
echo -e "\n" >> "formatted_born.dat"
awk '{if ($1 == 3) print $1, $2, $3, $4}' "$output_file" >> "formatted_born.dat"

echo "Processing complete; results saved to formatted_born.dat"

After extraction, simply sum the principal diagonal components for each element and take the average—this gives the effective Born charge Z* used in the THz force calculation.

3 MD simulation via Lammps

The MD protocol should follow the same setup used for standard BaTiO₃ calculations, 10 * 10 * 10 supercell is used for simulation, Temperature T is set to 300 K.

As with most Lammps simulations, you need to provide:

• Structure model as input file;
• force-field model as input file;
• appropriate MD parametersls: time step, Tdamp, Pdamp, T, etc.

One of the critical part is the application of the force acting on the atoms. That is, the electric field is converted into a force acting on the atoms — obtained by multiplying the Born effective charge Z* with the instantaneous electric field E(t) to give F = Z·E(t) —— and applied as an additional force on the ions. For example, when an electric field of Et=1MV/cm is applied, the force is F = 0.01 * (Z).

Another critical part is the THz electric field itself must be properly defined. We use a set of variables to fully characterize the THz electric field:

• Peak amplitude (V Å⁻¹ – must be converted according to the chosen LAMMPS units);
• THz carrier frequency (determines the time between two successive field maxima);
• Single-pulse centre time t₀ (s)；
• Pulse duration – Gaussian envelope standard deviation σ (s);
• Pulse interval – delay of 2 ps before the first pulse is fired (s);
• Total number of pulses to be delivered；
• Pulse counter – tracks how many pulses have already been applied.

These parameters together determine the shape and periodicity of the THz electric field which is corresponding to parameters in THz electric field is implemented as:

E(t)=Ef⋅exp[ −4 ln(2) (t−t₀)²/σ² ]⋅sin(2π⋅freq⋅t),

where

• Ef – peak amplitude (V Å⁻¹);
• t₀ – pulse center(ps);
• σ – full-width-at-half-maximum (FWHM) of the Gaussian envelope (ps);
• freq – carrier frequency (THz);
• t = (v_current_step - NSTEPS) * 1e-15 (elapsed simulation time);
• count is used to apply multiple pulses.

The Gaussian term confines the sinusoid to a single THz pulse, and the prefactor 4ln(2) guarantees that σ is the FWHM.

Additionally, we introduce v_current_step to keep track of the current timestep — this variable is updated every step and used to compute the instantaneous THz field and apply the corresponding force. The fore is defined as:

variable Bafx ${Z_Ba}*${Ef}*exp(-4*log(2)*(((v_current_step-750000)*1e-15-1000*${count}*1e-14-${t0})^2)/(2*(${sigma}^2)))*sin(2*PI*${freq}*((v_current_step-750000)*1e-15-1000*${count}*1e-14))

It is recommended to adopt the settings illustrated in the accompanying input.lammps template below. These variable parameters also can be adjusted as needed.

variable T_change index  300
variable Ef index  0.01 ## Peak amplitude (V Å⁻¹ – must be converted according to the chosen LAMMPS *units*)
variable freq index 1e12   # THz carrier frequency (determines the time between two successive field maxima)
variable t0 index 5e-12    # Single-pulse centre time t₀ (s)
variable sigma equal 1e-12    # Pulse duration – Gaussian envelope standard deviation σ (s)
variable interval equal 2e-12  # Pulse interval – delay of 2 ps before the first pulse is fired (s)
variable total_pulses index  1    # Total number of pulses** to be delivered
variable count equal 0          # Pulse counter – tracks how many pulses have already been applied



shell mkdir ${t0}
shell cd ${t0}
shell cp ../input.lammps ./
shell cp ../input/*.pb ./
shell cp ../input/*.restart ./
shell cp ../input/input.lammps ./
log log.${T_change}

variable        TEMP        equal ${T_change} 
variable        NSTEPS          equal 750000
variable        THERMO_FREQ     equal 10
variable        DUMP_FREQ       equal 10
variable        PRES            equal 0.000000
variable        TAU_T           equal 0.100000
variable        TAU_P           equal 12
variable    time_step   equal 0.001

units           metal
boundary        p p p
atom_style      atomic

neighbor        1.0 bin

#read_data       200.data
read_restart    BTO_101010.restart
change_box      all triclinic

mass            1 137.327000
mass            2 207.200000
mass            3 40.078000
mass            4 87.620000
mass            5 208.980400
mass            6 39.098300
mass            7 22.989700
mass            8 178.490000
mass            9 47.867000
mass            10 91.224000
mass            11 92.906000
mass            12 23.305000
mass            13 114.818000
mass            14 65.409000
mass            15 15.999400

plugin load libdeepmd_lmp.so
pair_style      deepmd ./graph.pb
pair_coeff      * *

#### apply Electric Field ####
variable Z_Ba equal 2.77
variable Z_Ti equal 6.26
variable Z_O equal -3.01


group  Ba type 1
group  Ti type 9 
group  O  type 15



label pulse_loop
timestep 0.001
thermo          ${THERMO_FREQ}



dump           ${count}_1  all custom ${DUMP_FREQ} ${count}traj.lammpstrj id type x y z 
dump_modify    ${count}_1   sort id
#velocity        all create ${TEMP} 114514
variable step_time  equal ${count}*${time_step} 
variable current_step equal step
variable count_interval equal ${count}*${interval} 

#variable center equal ${step_time}-${count}*${interval}-${t0}
#variable shifted_step equal ${step_time}-${count}*${interval}
variable inercos equal 2*PI*${freq}*((step-750000)*1e-15*10-${count_interval})

variable VALUEZ  equal 0
variable VALUEY  equal 0
#variable VALUEZ equal ${Ef}*exp(-4*log(2)*(((step-750000)*10*1e-15-${count_interval}/10-${t0})^2)/(2*(${sigma}^2)))*sin(2*PI*${freq}*((step-750000)*10*1e-15-${count_interval}/10))
variable VALUEX equal ${Ef}*exp(-4*log(2)*(((v_current_step-750000)*1e-15-1000*${count}*1e-14-${t0})^2)/(2*(${sigma}^2)))*sin(2*PI*${freq}*((v_current_step-750000)*1e-15-1000*${count}*1e-14))
#variable VALUEZ equal ${Ef}*exp(-4*log(2)*(((step-750000)*10*1e-15-${t0})^2)/(2*(${sigma}^2)))*sin(2*PI*${freq}*((step-750000)*10*1e-15))


#variable V_exp equal exp(-4*log(2)*(${center}^2)/(2*(${sigma}^2)))

#thermo_style    custom step temp v_step_time  v_count etotal press vol lx ly lz v_VALUEX v_VALUEY v_VALUEZ


variable Bafx equal ${Z_Ba}*${Ef}*exp(-4*log(2)*(((v_current_step-750000)*1e-15-1000*${count}*1e-14-${t0})^2)/(2*(${sigma}^2)))*sin(2*PI*${freq}*((v_current_step-750000)*1e-15-1000*${count}*1e-14))
variable Bafy equal ${Z_Ba}*${VALUEY}
variable Bafz equal ${Z_Ba}*${VALUEZ}
variable Tifx equal ${Z_Ti}*${Ef}*exp(-4*log(2)*(((v_current_step-750000)*1e-15-1000*${count}*1e-14-${t0})^2)/(2*(${sigma}^2)))*sin(2*PI*${freq}*((v_current_step-750000)*1e-15-1000*${count}*1e-14))
variable Tify equal ${Z_Ti}*${VALUEY}
variable Tifz equal ${Z_Ti}*${VALUEZ}
variable Ofx equal ${Z_O}*${Ef}*exp(-4*log(2)*(((v_current_step-750000)*1e-15-1000*${count}*1e-14-${t0})^2)/(2*(${sigma}^2)))*sin(2*PI*${freq}*((v_current_step-750000)*1e-15-1000*${count}*1e-14))
variable Ofy equal ${Z_O}*${VALUEY}
variable Ofz equal ${Z_O}*${VALUEZ}

fix fBa Ba addforce v_Bafx v_Bafy v_Bafz
fix fTi Ti addforce v_Tifx v_Tify v_Tifz
fix fO  O  addforce v_Ofx v_Ofy v_Ofz
print "Bafx = ${Tifx}"
print "Bafx = ${Tify}"
print "Bafx = ${Tifz}"
print "Current step = ${current_step}"

compute force_${count}Ti Ti property/atom fx fy fz
dump Ti_2_${count} Ti custom 100 dump.force_Ti id  type x y z c_force_${count}Ti[1] c_force_${count}Ti[2] c_force_${count}Ti[3]                   # You can choose to output the force in the corresponding direction.
dump_modify   Ti_2_${count}   sort id
thermo_style    custom step temp v_step_time  v_count etotal press vol lx ly lz v_Bafx v_Bafy v_Bafz


#fix             monitor all print 100 "TIME=${count} EZ=${VALUEZ} Step_time=${step_time}  Cen=${center} Shift=${shifted_step} V_exp=${V_exp} " file ef_profile.dat   
fix             100 all npt temp ${TEMP} ${TEMP} ${TAU_T} tri ${PRES} ${PRES} ${TAU_P}
run         10000                       # Run for 10 ps (corresponding to interval)）
unfix       fBa                         # remove the electric field
unfix       fTi                         # remove the electric field
unfix       fO                          # remove the electric field
unfix   100
variable    count equal ${count}+1      # Increment the counter by one
print "Current count: ${count}"
if          "${count} < ${total_pulses}" then "jump SELF pulse_loop"
write_data      efield${TEMP}.data 
write_restart   efield${TEMP}.restart
shell cd ..
clear 
next t0
jump input.lammps

After this run you will obtain

• log.300 – LAMMPS log (thermo output)
• 300trj.lammpstrj – corresponding trajectory file

The supercell volume and temperature written in log.300, as well as the position information stored in the .lammpstrj trajectory file, and polarization of different steps can be extracted automatically with the supplied parsing script from Li Denan's markdown.(https://github.com/MoseyQAQ/FerroDispCalc/tree/cpp).

#Px Py Pz
Px1 Py1 Pz2  #step 1
Px2 Py2 Pz2  #step 2
...          
Pxn Pyn Pzn  #step n

Once we have saved the polarization data file (polar.dat), we can immediately use it for plotting.