GCMC with molecules

111291 · July 28, 2016, 3:48pm

Dear all,

I am trying to run gcmc with a molecule, ethane. I already tried with monatomic gas such as CH4 using molecule command and worked well.

However, I changed it into ethane (C2H6) and it did not work properly. This is first time to try gcmc with the molecules and I might miss something.

The results showed no change in the number of molecules, only reduction converging to zero, or lost atoms.

I am attaching three files, a input file and two data files.

I would appreciate any help.

Best regards,

Seunghwan Baek

Input file

variable t equal 353

variable p equal 137.90

variable seed equal 31798.0000

variable phi equal 0.4999

variable x equal 50.0

units metal

atom_style full

boundary p p p

pair_style lj/cut 14.0

bond_style harmonic

read_data C2_bulk_10.txt

molecule ethane C2H6.txt

special_bonds lj 1.0 0.0 0.0 angle no dihedral no

group gas type 1

region cv block 0.0 $x 0.0 $x 0.0 $x units box side in

pair_modify mix arithmetic

pair_modify tail yes

neighbor 1.5 bin

neigh_modify every 1 delay 0

velocity gas create t {seed} mom yes rot yes dist gaussian

minimize 1.0e-8 1.0e-8 100000 1000000

compute allperatom gas stress/atom NULL

compute allp gas reduce sum c_allperatom[1] c_allperatom[2] c_allperatom[3]

variable allpress equal -(c_allp[1]+c_allp[2]+c_allp[3])/(3*vol)

compute mdtemp gas temp

compute_modify mdtemp dynamic yes

fix sh gas shake 1.0e-8 200 0 b 1 t 1 m 15.035 mol ethane

fix GC gas gcmc 100 200 500 0 ${seed} $t -4.607 0.05 mol ethane shake sh region cv pressure p fugacity_coeff {phi} group gas full_energy

variable nm equal count(gas)/2

thermo_style custom step vol v_nm c_mdtemp v_allpress etotal

thermo_modify format float %24.8g

thermo 1

timestep 0.0001

run 1000

Data file - initial configuration

LAMMPS data file. C2_test.

20 atoms

10 bonds

0 angles

0 dihedrals

0 impropers

1 atom types

1 bond types

0 angle types

0 dihedral types

0 improper types

0.0000 50 xlo xhi

0.0000 50 ylo yhi

0.0000 50 zlo zhi

Masses

1 15.035000 # CH3(sp3)

Pair Coeffs

1 0.00844 3.7500 # CH3(sp3)

Bond Coeffs

1 13.44 1.54 # CH3-CH3

Atoms

1 1 1 0 2.22600000000000 1.00100000000000 1.00100000000000

2 1 1 0 1.00100000000000 1.90600000000000 1.26200000000000

3 2 1 0 12.2280000000000 1.00100000000000 1.00100000000000

4 2 1 0 11.0030000000000 1.90600000000000 1.26200000000000

5 3 1 0 17.2290000000000 1.00100000000000 1.00100000000000

6 3 1 0 16.0040000000000 1.90600000000000 1.26200000000000

7 4 1 0 22.2300000000000 1.00100000000000 1.00100000000000

8 4 1 0 21.0050000000000 1.90600000000000 1.26200000000000

9 5 1 0 27.2310000000000 1.00100000000000 1.00100000000000

10 5 1 0 26.0060000000000 1.90600000000000 1.26200000000000

11 6 1 0 32.2320000000000 1.00100000000000 1.00100000000000

12 6 1 0 31.0070000000000 1.90600000000000 1.26200000000000

13 7 1 0 37.2330000000000 1.00100000000000 1.00100000000000

14 7 1 0 36.0080000000000 1.90600000000000 1.26200000000000

15 8 1 0 42.2340000000000 1.00100000000000 1.00100000000000

16 8 1 0 41.0090000000000 1.90600000000000 1.26200000000000

17 9 1 0 2.22600000000000 6.80200000000000 1.00100000000000

18 9 1 0 1.00100000000000 7.70700000000000 1.26200000000000

19 10 1 0 12.2280000000000 6.80200000000000 1.00100000000000

20 10 1 0 11.0030000000000 7.70700000000000 1.26200000000000

Bonds

1 1 1 2

2 1 3 4

3 1 5 6

4 1 7 8

5 1 9 10

6 1 11 12

7 1 13 14

8 1 15 16

9 1 17 18

10 1 19 20

Molecule information

LAMMPS C2H6 molecule file

2 atoms

1 bonds

0 angles

0 dihedrals

0 impropers

30.070 mass

0.0 0.0 0.0 com

0.0 0.0 0.0 0.0 0.0 0.0 inertia

Coords

1 2.22600000000000 1.00100000000000 1.00100000000000

2 1.00100000000000 1.90600000000000 1.26200000000000

Types

1 1

2 1

Masses

1 15.035

2 15.035

Bonds

1 1 1 2

Special Bond Counts

1 1 0 0

2 1 0 0

Special Bonds

1 2

2 1

Shake Flags

1 2

2 2

Shake Atoms

1 1 2

2 1 2

Shake Bond Types

1 1

2 1

C2_bulk_10.txt (1.61 KB)

C2H6.txt (485 Bytes)

main.in (1.41 KB)

athomps · July 28, 2016, 11:21pm

There are many reasons why your simulation is not behaving as expected:

"The results showed no change in the number of molecules, only reduction converging to zero, or lost atoms. "

Maybe you set the chemical potential too low.
Maybe your molecule template is very high energy and/or high forces

The key is to try to simplify your problem to something that works, and then incrementally add things to get to where you want to go. The first step would be to test your molecule template using create_atoms command. See if you can run NVE MD with a few molecules.

111291 · July 30, 2016, 12:45am

Thank you, Aidan,

I tried to change my chemical potential (-50 ~ -0.10), which gave the same result (reduction converging to zero).

I think that my chemical potential comes from right calculation and also used pressure command. I used to run NVE, NVT or NPT under this pressure.

In addition, I tried to run simple NVE with the data file, not with “create_atoms” command, having 10 molecules and it ran.

Step Volume nm mdtemp allpress TotEng

54 125000 10 428.10638 -150984.72 49.358547

100000 125000 10 566.1166 5.2240923 14604.103

200000 125000 10 572.11607 2.9483276 14604.103

300000 125000 10 566.87224 4.7355664 14604.103

400000 125000 10 565.78646 5.9968002 14604.103

500000 125000 10 565.77721 5.3452818 14604.103

600000 125000 10 567.19402 5.3180585 14604.103

700000 125000 10 575.3314 5.5850096 14604.103

800000 125000 10 566.94418 4.7988622 14604.103

900000 125000 10 565.98081 5.0106637 14604.103

1000000 125000 10 566.84835 4.6082831 14604.103

1000054 125000 10 567.01378 4.6020261 14604.103

Loop time of 15.8189 on 20 procs (20 MPI x 1 OpenMP) for 1000000 steps with 20 atoms

What else can I try?

Seunghwan Baek

athomps · July 30, 2016, 2:27am

Do this:

“The first step would be to test your molecule template using create_atoms command. See if you can run NVE MD with a few molecules.”

111291 · August 4, 2016, 1:34pm

Dear Aidan,

I have run NVE with create_atoms command and I found what I missed.

Another question is that when I run fix gcmc (Before GCMC, fix NVT was run for initial equilibrium), the temperature and pressure does not change into set points of them.

what causes this?

Seunghwan

athomps · August 4, 2016, 7:07pm

Is this not clear?

“Inserted atoms and molecules are assigned random velocities based on the specified temperature T. Because the relative velocity of all atoms in the molecule is zero, this may result in inserted molecules that are systematically too cold. In addition, the intramolecular potential energy of the inserted molecule may cause the kinetic energy of the molecule to quickly increase or decrease after insertion. The tfac_insert keyword allows the user to counteract these effects by changing the temperature used to assign velocities to inserted atoms and molecules by a constant factor. For a particular application, some experimentation may be required to find a value of tfac_insert that results in inserted molecules that equilibrate quickly to the correct temperature.”

111291 · September 16, 2016, 9:36pm

Thank you, Aidan

With the keyword, tfac_insert, the temperature fluctuated a little bit around the set point. However, the pressure is still far away from what I want. I have tried to change the values of the keywords, pressure and fugacity_coeff but it seemed not to work out.

The number of molecules here changed forward around 1100. When I calculated with the ideal gas law, it looks similar. This means that C2 comes in properly and the pressure calculation is wrong. This is bulk state simulation and the compute command and thermo_style(press) give the same values.

Do you have any suggestion?

compute allperatom gas stress/atom NULL

compute allp gas reduce sum c_allperatom[1] c_allperatom[2] c_allperatom[3]

variable allpress equal -(c_allp[1]+c_allp[2]+c_allp[3])/(3*vol)

compute mdtemp gas temp

compute_modify mdtemp dynamic yes

fix MC gas gcmc 500 700 700 0 31798 353 -2.819 0.05 mol ethane shake const_bond region substrate pressure 137.9 fugacity_coeff 0.4999 group gas full_energy tfac_insert 1.80

WARNING: Molecule attributes do not match system attributes (…/molecule.cpp:1350)

compute_modify mdtemp dynamic yes

variable num_C2 equal count(gas)/2

thermo 10000

thermo_style custom step v_num_C2 c_mdtemp v_allpress press

WARNING: New thermo_style command, previous thermo_modify settings will be lost (…/output.cpp:680)

thermo_modify format float %24.16g

neighbor 1.5 bin

neigh_modify every 1 delay 0

run 100000

0 atoms in group FixGCMC:gcmc_exclusion_group:MC

0 atoms in group FixGCMC:rotation_gas_atoms:MC

Neighbor list info …

1 neighbor list requests

update every 1 steps, delay 0 steps, check yes

max neighbors/atom: 2000, page size: 100000

master list distance cutoff = 15.5

ghost atom cutoff = 15.5

binsize = 7.75 -> bins = 7 7 7

Memory usage per processor = 9.43478 Mbytes

Step num_C2 mdtemp allpress Press

0 1 339.0829455701206 -2.07269953325158e-06 -8.944009075565275e-05

10000 823 376.7174762557219 -21617.23298558185 -21617.23307294924

20000 946 376.0856606501425 -24709.0862059211 -24709.0862932885

30000 1007 375.2185382916839 -26016.15480931733 -26016.15489668473

40000 1056 376.2512756013458 -27315.58032180904 -27315.58040917644

50000 1084 376.3326795500135 -27814.9245516896 -27814.92463905699

60000 1093 378.4262384389432 -28014.11331053155 -28014.11339789894

70000 1099 379.6100605253111 -28289.45024187829 -28289.45032924568

80000 1109 381.4981055796842 -28460.8075657948 -28460.80765316218

90000 1110 382.508423202094 -28091.92873479091 -28091.9288221583

100000 1123 384.1637239308351 -28660.23640896524 -28660.23649633263

Loop time of 2562.39 on 20 procs for 100000 steps with 2246 atoms

111291 · September 18, 2016, 4:28pm

Dear all,

I am trying to run gcmc with a molecule, ethane. The pressure is still far away from what I want and the number of ethane is too high. I have tried to change the values of the keywords, pressure and fugacity_coeff but it seemed not to work out. The density is higher than the reference density in 2000 psi but pressure shows negative.

I need your advice!