Dear all,
I am using the potential by M. Sprik et al. ( J.phys. Chem. 96 2027 (1992) ) to simulate a C60 molecular crystal. In this framework the C60 molecules are treated as rigid bodies made of 90 interaction centers (60 carbon atoms + 30 sites in the middle of each double bond) the sites interact both via a short-range Lennard-Jones potential and via a long-range Coulomb term. Unfortunately even at zero temperature I cannot recover both the lattice spacing and the correct total energy shown in the Sprik paper. I am almost sure the error as to do with the Ewald summation.
Here is the input file, the unit cell is a simple cubic one, with lattice parameter 14.04 and a basis of 4 neutral molecules (360 atoms, the complete structure file is attached):
units real
atom_style full
kspace_style ewald 1.0e-4
boundary p p p
read_data initpos.dat
mass 1 12.011
mass 2 0.0001
velocity all create 0.00000 87287
pair_style lj/cut/coul/long 14.04078 14.04078
pair_coeff 1 1 0.0298 3.4
pair_coeff 1 2 0.0298 3.5
pair_coeff 2 2 0.0298 3.6
neighbor 0.3 bin
neigh_modify every 20 delay 0 check no one 50000 page 5000000
neigh_modify exclude molecule all
fix 1 all rigid molecule
thermo 1
thermo_style custom step temp ke pe etotal
run 10
- changing the Ewald sum accuracy from 10e-4 to 10e-12 the total energy ranges from -380 Kcal/mol to -1500 Kcal/mol, according to the Sprik paper it should be around -180 Kcal/mol
- keeping fixed the accuracy at 10-10 and changing the cut-off radius, the total energy ranges from -759 Kcal/mol to -610 Kcal/mol when the cut-off is changed from 14 A to 40 A
To my understanding, in the standard Ewald theory, the total energy should converge to a stable value once the cut-off and the accuracy are chosen large enough. Can anybody help me? Am I missing some important settings?
Thanks in advance, best regards.
Andrea
initpos.dat (15.5 KB)