Inconsistency between Lammps, GOMC, and and Cassandra at High Densities

Hi everyone,

1. I have performed three sets of NVT simulations in LAMMPS, GOMC, and Cassandra. My system contains 300 molecules of n-butane at 510K and TraPPE-UA force field was used for all simulations. The densities range from 0.1 gcc to 0.65 gcc.

At low densities, pressures seem to match perfectly, but at high densities, there is a significant difference between pressures computed by these three codes. Cassandra is further away from LAMMPS than GOMC. I have been communicating with GOMC people and I've made sure that my simulations are correct, and we decided that the inconsistency between GOMC and LAMMPS might be due to the fact that GOMC still lacks regrowth move. Since I'm new to Cassandra, I would like to make sure that I have done the simulations correctly, and there is no technical mistakes in my input files, even though I have used the sample simulations to make my input files. Attached is one of my simulations at rho=0.662 gcc and T=510K, along with a plot of three isotherms simulated by Lammps, GOMC and Cassandra.


I also have these a few other related questions that I couldn't find an answer for in the documentation:

2. Does Energy_LJ include both intermolecular and intramolecular LJ interactions, or only intermolecular?

3. Does Cassandra neglect 1-4 bonded LJ interactions by default? I did not specify anything for Intra_Scale section and Cassandra didn't compain, so I'm assuming it must have a default value.




Thank you,

-Mostafa Razavi

Hi Ryan,

Thank you for your helpful response.

1- Following are the pressures in atm that I computed. You can also see the attached plot for a comparison:

LAMMPS: 2671 +/- 24
Cassandra: 2742 +/- 214
GOMC: 2695 +/- 80

Your results:
Gromacs: 2755 +/- 304
Cassandra: 2727 +/- 221

I run Cassandra simulations for 90000 sweeps and started averaging after 30000 sweeps.

There is a big difference between my LAMMPS and your Gromcs results. I used everything according on TraPPE except I have harmonic bond (E=K*(r-r0)^2) and K=191.75 kcal/mol/A^2. What was you Kbond value? Can you send me your Gromacs files so I can compare with mine?

3- I actually did try the following

# Intra_Scaling
0. 0. 0. 1.
0. 0. 0. 1.

However, the results were identical to when I used default values. Do you know why? I think I will double check though.

Thanks,

-Mostafa

OK now it makes more sense. Your Cassandra results (energy-lj and pressure) matches GOMC perfectly. My Cassandra results got closer when I excluded 1-4 interactions (I had missed the underscore in Intra_Scaling, so Cassandra did not see it and ironically did not give me error) and chose "equilibration 1100" for my Run_Type 1-4, but is still a little off. I wonder what I am missing in my input files. Could you send me your Cassandra simulation files at 510K and 0.662 gcc, or could you please compare my attached most recent Cassandra files with yours?

FYI: I have a hard time uploading files in here. I have to zip everything and it shouldn't exceed 0.1MB. That's why I couldn't attach the log file. I can email you though.

Regarding the uncertainties, I took the average of the values reported in nvt.out.prp after 30000 sweeps and calculated the standard deviations for them. Same goes for GOMC. In LAMMPS I took the average of instantaneous values.

I know it's a long shot, but is there any chance you could repeat your Gromacs simulation using harmonic bond (E=K*(r-r0)^2) and K=191.75 kcal/mol/A^2? We are studying the performance of MC vs MD in vapor pressure calculations and your help would be very much appreciated.


Thank you

-Mostafa