Unsatisfying isotropic phase with Gay-Berne NPT

Dear all,

I am trying to reproduce the phase diagram of discotic liquid crystals (oblate ellipsoids) using Gay-Berne potential, but failed at the beginning to generate the isotropic phase.

[1] Phys. Rev. E 67, 041703 (2003) - Influence of shape and energy anisotropies on the phase diagram of discotic molecules (MC)
[2] Cookie Absent (MD,LAMMPS)
For shape anisotropic ratio 0.2 (1:5) and energy ratio 0.1 (1:10), isotropic state at T=12, P=100. [1][2]
For shape anisotropic ratio 0.345 and energy ratio 0.2, isotropic at T=4, P=50. [1]

My run directly go to a more orientationally ordered phase, also the internal energy/enthalpy is lower than expected. The internal energy per particle is around 90 in Busselez 2014, but mine is already negative.

I use npt/asphere to control the system under reduced unit following [2] the MD study using LAMMPS, all units and timestep are set the same.
The energy unit is directly 1 at the crossed configuration [1],also checked with two single particles using the same pair style.

Since the potential setting should be correct, I cannot tell if the NPT setting is wrong.
Please tell me which part I should fix to get the isotropic phase.


Quick script for T=12, P=100 case:

#repeat aspect ratio 0.2 with Busselez 2014
units		lj
atom_style	ellipsoid
atom_modify	first big
dimension	3

# create big ellipsoidal particles

lattice		sc 0.1
region		box block 0 16 0 16 0 16
create_box	2 box
create_atoms	1 region box

set		type 1 mass 1.0
set		type 1 shape 1.0 1.0 0.2
group		big type 1
set		group big quat/random 29890

velocity	big create 24.0 87280 mom yes loop all

# equilibrate big particles
variable        ratio equal 0.2*0.2
variable        xxx equal (${ratio}-1)/(${ratio}+1)
print           ${xxx}
variable        e0 equal (1-${xxx}*${xxx})
print           ${e0}

pair_style      gayberne/omp 1.0 2.0 1.0 2.0
pair_coeff      * * 1.0 1.0 1.0 1.0 10.0 0 0 0

neighbor	1 bin
neigh_modify	delay 0 every 1 check yes one 4000

# compute temperature and pressure, setting NPT/ASPHERE

compute         temp_trans all temp

variable        dof equal count(big)
compute         temp_rottrans all temp/asphere dof all
compute_modify  temp_rottrans extra/dof ${dof}

compute         pres_trans all pressure temp_trans
compute         pres_ke all pressure temp_trans ke
compute         pres_virial all pressure temp_trans virial

compute         p_e all pe
compute         ke_t all ke
compute         ke_r all erotate/asphere

variable        erot equal c_ke_r/4096

variable        volume equal vol

compute         orient all property/atom quati quatj quatk quatw
compute         diameter all property/atom shapex shapey shapez 

# output in time average and log file

fix		1 big npt/asphere/omp temp 12.0 12.0 0.025 iso 100.0 100.0 0.025
fix_modify      1 temp temp_rottrans
fix_modify      1 press pres_trans
#compute_modify  1_temp extra/dof ${dof}
fix             3 all ave/time 1 100 100 v_volume c_temp_rottrans c_pres_trans c_p_e c_ke_t c_ke_r file temp12_0.avg

dump		1 all custom 1000 temp12_0.dump id type x y z &
                   c_orient[1] c_orient[2] c_orient[3] c_orient[4] &
                   c_diameter[1] c_diameter[2] c_diameter[3]

thermo_style	custom step c_temp_rottrans press etotal pe ke v_erot enthalpy 
thermo		1000


timestep        0.00025
run		500000

write_restart   temp12_0.restart
undump          1
unfix           1
unfix           3

Why is the reason for not including the kinetic energy in the calculation of the pressure? This is a tricky point in MD simulations of 3D objects, and I came up with a similar solution to yours for thermostatting and barostatting. My understanding is that the pressure for rigid bodies should include the translational component of the temperature (as you have done) and the full virial, including the kinetic energy, i.e.

compute  press_transl all pressure temp_trans
fix		1 big npt/asphere/omp temp 12.0 12.0 0.025 iso 100.0 100.0 0.025
fix_modify      1 temp temp_rottrans
fix_modify      1 press press_transl

Let us know if that works for you.

Thank you for your reply.
I actually set the thermo/barostatting like this after I read another discussion topic between you and sb. on NPT/ASPHERE. Especially the temperature part should not include the rotational part of energy when controlling pressure.
for the code of compute pressure It use

scalar = (temperature->dof * boltz * temperature->scalar +
                virial[0] + virial[1] + virial[2]) / 3.0 * inv_volume * nktv2p

If the particle is asphere, this temperature->dof term will not be 3N.
this will not be a problem when the temperature is based on translational.
But another point is that the translational temperature will only be equal to rotational temperature at equilibrium. if the system is not at equilibrium, the value of temp and temp/asphere will not be the same, which may be a problem?

In my script I use:

compute         temp_trans all temp
compute         pres_trans all pressure temp_trans
fix		        1 big npt/asphere/omp temp 12.0 12.0 0.025 iso 100.0 100.0 0.025
fix_modify      1 temp temp_rottrans
fix_modify      1 press pres_trans

this is the same as in your example, only the name of the compute is different, is this right?
This works for my pressure control.

When using a different setting, the npt/asphere setting in some examples

fix		        1 big npt/asphere/omp temp 12.0 12.0 0.025 iso 100.0 100.0 0.025
compute_modify  1_temp extra/dof ${dof}

It cannot (may not?) reach the target pressure and sth else went wired when I monitor the component pxx, pyy and pzz in a different simulation . So I think the translational temperature to pressure one may be better.

Yes, my bad! I overlooked your script regarding which pressure was fed to the barostat.

I don’t think this is the problem, as the thermostat will rapidly populate both degrees of freedom. Let me have a look at your simulation again and get back to you.

I think the problem is in the ambiguity in which the GB parameters are defined. In the papers you cited, the authors give only relative values such as k'=0.1, defined as k'=\varepsilon_\text{ee}/\varepsilon_\text{ff}. In your specific implementation, you set \varepsilon_\text{ff}=10, i.e.

pair_coeff      * * 1.0 1.0 1.0 1.0 10.0 0 0 0

A more formal way to input these parameters being:

pair_coeff      * * 10.0 1.0 0.1 0.1 1.0 0.1 0.1 1.0

The difference in the phase diagram arises from the depth of the potential energy well. If you set the well depth to 1.0, you get the isotropic phase at T=12* and P=100*, as you can see in the render. However, you better check a few more points in the phase diagram of Fig 10/d of reference [1].

Sorry for the off-topic, but what do you use for rendering?

The image was rendered with Intel OSPRay, as implemented in Ovito. I have added the deep of field with an aperture of 0.4 to slightly defocus the scene.

1 Like

Hi hothello,
Thank you for your help!
I think I did not understand the GB parameter well.
I also consult the author of reference [1], she suggest beside the ones you point out, the effective radius parameter should set to be 0.2, which is the thickness of disc particle:

pair_coeff      * * 10.0 0.2 0.1 0.1 1.0 0.1 0.1 1.0

I get the isotropic phase I want with the corrected pair coefficient setting.

1 Like