I would like to simulate a dipolar Gay-Berne potential to study the liquid crystalline properties of the system. Can you guide me on how to modify the GB potential to incorporate dipole interactions? I’m currently using this script for the system, but the temperature is not being controlled properly during the simulation. Could you help me identify what went wrong in my simulation process?
variables
------------------------------------------------------------------------------
----------- CONSTANT VARIABLES -----------------------------------------------
------------------------------------------------------------------------------
variable PI equal (3.1415926535897932384626433832795)
variable SQRT32 equal (0.86602540378443864676)
variable SQRT3 equal (1.73205080756887729353)
------------------------------------------------------------------------------
----------- SYSTEM VARIABLES -------------------------------------------------
------------------------------------------------------------------------------
variable T equal 2.00
variable P equal 1.20
variable rho equal 0.05
variable nm equal (5)
variable nm2 equal (v_nm/2)
variable as equal (sqrt(1.0/v_rho))
variable Lx equal (v_nm*v_as)
variable Ly equal (v_Lx)
variable Lz equal (v_Lx) # Added for 3D system
variable N equal (v_nmtot)
variable stime equal 400000
variable rtime equal 2000000
------------------------------------------------------------------------------
----------- Gay-Berne ellipsoids with dipole moments in a 3D box --------------
------------------------------------------------------------------------------
units lj
atom_style hybrid ellipsoid dipole # Modified for dipole moments
#atom_style hybrid molecular ellipsoid dipole
dimension 3 # Changed to 3 for 3D system
#read_restart restart.gb.no.subs.npt.sat.P.${P}.2
lattice sc {rho} # Changed to simple cubic for 3D system
region box block 0 {nm} 0 {nm} 0 {nm} # Modified for 3D
create_box 1 box
create_atoms 1 box
set type 1 mass 1.0
set type 1 shape 3 1 1
set group all quat/random 123456782
Set random dipole moments for each particle
set type 1 dipole/random 12345 1.1 # Assign random dipole moments
variable gb_ps_gamma equal 1.0
variable gb_ps_upsilon equal 1.0
variable gb_ps_mu equal 2.0
variable gb_ps_cutoff equal 5.0
variable dipole_cutoff equal 1.0
variable gb_pc_epsilon equal 1.0
variable gb_pc_sigma equal 1.0
variable gb_pc_epsilon_i_a equal 1.0
variable gb_pc_epsilon_i_b equal 1.0
variable gb_pc_epsilon_i_c equal 0.2
variable gb_pc_epsilon_j_a equal 0.0
variable gb_pc_epsilon_j_b equal 0.0
variable gb_pc_epsilon_j_c equal 0.0
variable gb_pc_cutoff equal 5.0
pair_style hybrid/overlay gayberne 1.0 2.0 1.0 4.0 lj/cut/dipole/long 4.0
#pair_style hybrid/overlay lj/cut/dipole/long 3.0 gayberne 1.0 2.0 1.0 5.0
pair_coeff * * gayberne 1.0 1.0 1.0 1.0 0.2 1.0 1.0 0.2
pair_coeff * * lj/cut/dipole/long 0.0 0.0
kspace_style ewald/dipole 0.00001
------------------------------------------------------------------------------
----------- Prepare the Ensemble ---------------------------------------------
------------------------------------------------------------------------------
compute rot all temp/asphere
group uniaxial type 1
variable dof equal count(uniaxial)
compute_modify rot extra/dof ${dof}
compute pr all pressure rot
velocity all create ${T} 987654329 mom yes rot yes dist gaussian
neighbor 0.8 bin
#neigh_modify skin 3.0
neigh_modify delay 0
timestep 0.002
Modify thermo_style to output dipole properties
thermo_style custom step c_rot temp epair etotal c_pr press vol density
thermo 1000
compute q all property/atom quatw quati quatj quatk
#min_modify line quadratic
#min_style cg
#minimize 0.01 1.0e-05 10000 100000
#minimize 0.00 1.0e-10 10000 100000
------------------------------------------------------------------------------
---------------------- NVT Ensemble ------------------------------------------
------------------------------------------------------------------------------
#undump 1
dump 1 all custom 10000 dump.gb.no.subs.nvt.sat.P.${P} id x y z c_q[1] c_q[4] fx fy fz
#dump_modify 1 format line “%d %20.15e %20.15e %20.15e %20.15e %20.15e %20.15e”
#unfix 2
fix 2 all nvt/asphere temp 1.40 1.40 0.1
#fix 2 all nvt/asphere temp {T} {T} 0.1 mtk no pchain 0 tchain 1
compute_modify 2_temp extra/dof ${dof}
restart 100000 restart.gb.no.subs.nvt.sat.P.{P}.1 restart.gb.no.subs.nvt.sat.P.{P}.2
run ${stime}
------------------------------------------------------------------------------
---------------------- Data collection run -----------------------------------
------------------------------------------------------------------------------
------------------------------------------------------------------------------
------------------------------------------------------------------------------
compute RDF all rdf 500 cutoff 4.0
fix RDF_OUTPUT all ave/time 25 100 5000 c_RDF[*] file rdf_0.1_lj.out mode vector
run 5000