[gmx-users] Problem in CaF2 with Buckinghan Force Field

[Tandia, Adama]

Dear All,
I'm running a MD NPT for CaF2 with Buckingham force field. The force field parameters are known to be very good. I got them from literature and many people had used them in the past. When I run the MD for let us say 50 ps, the density of the system goes to zero (0.02...)if the run does not crash, and it stays so no matter how long you run the MD after. There are may be some thing that I'm missing. I'm using the windows version 3.1.2 of GROMACS.
Has anybody ever succeeded on running MD with Buckingham for condensed phase materials? I search the list and found that earlier some people had had trouble to reproduce data with Silica ...
I enclose the input files (.mdp, .top).
I tried many different combinations of variables in the Electrostatic part, but nothing really changed.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
[ defaults ]
; non-bond type  comb-rule   
2 		     1	   
[ atomtypes ]
;name    	mass      	charge   	ptype   	A	rho	C6		
F		18.998	0.0		A	      0.0	0.0	0.0
Fn		18.998	0.0		A		0.0	0.0	0.0
Ca		40.078	0.0		A		0.0	0.0	0.0
		
[ nonbond_params ]
;i   j     func      	a                	b       	c6
F    Ca    2		65061.1841		29.76190476	0.0
F    F     2		174448.496		34.12969283	0.010526732								
Ca   Ca    2		0.0			 0.0010001	0.0

Fn   Ca    2		65061.1841		29.76190476	0.0
Fn   Fn    2		174448.496		34.12969283	0.010526732		
F    Fn    2		174448.496		34.12969283	0.010526732

[ moleculetype ]
; name nrexcl
F     0
[ atoms ]
; nr  type resnr residue 	atom  cgnr   charge     mass
1     F    1     F      	F    	1      -1.0    	18.998 

[ moleculetype ]
; name nrexcl
Fn     0
[ atoms ]
; nr  type resnr residue 	atom  cgnr   charge     mass
1     Fn   1     Fn      	Fn   	1      0.0    	18.998 

[ moleculetype ]
; name nrexcl
Ca	  0
[ atoms ]
; nr  type resnr residue 	atom  cgnr   charge     mass
1     Ca   1     Ca           Ca    1      2.0        40.078
				
[ system ]
 Bulk

[ molecules ]
F 512
Ca 256

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; VARIOUS PREPROCESSING OPTIONS = 
title                    = Yo
cpp                      = C:\gromacs\bin\cpp
; RUN CONTROL PARAMETERS = 
integrator               = md
dt                       = 0.00005
nsteps                   = 100000
; number of steps for center of mass motion removal = 
nstcomm                  = 1

; OUTPUT CONTROL OPTIONS = 
; Output frequency for coords (x), velocities (v) and forces (f) = 
nstxout                  = 1000
nstvout                  = 1000
nstfout                  = 1000
; Output frequency for energies to log file and energy file = 
nstlog                   = 5000
nstenergy                = 1000
; Output frequency and precision for xtc file = 
nstxtcout                = 5000
xtc_precision            = 5000
; This selects the subset of atoms for the xtc file. You can = 
; select multiple groups. By default all atoms will be written. = 
;xtc-grps                 = 
; Selection of energy groups = 
energygrps               = F Ca 

; NEIGHBORSEARCHING PARAMETERS = 
; nblist update frequency = 
nstlist                  = 10
; ns algorithm (simple or grid) = 
ns_type                  = simple
; Periodic boundary conditions: xyz or none = 
pbc                      = xyz
; nblist cut-off         = 
rlist                    = 0.9
; OPTIONS FOR ELECTROSTATICS AND VDW = 
; Method for doing electrostatics = 
coulombtype              = Cut-off
rcoulomb                 = 0.9
rvdw                     = 0.9
fourierspacing           = 0.12
pme_order                = 4
ewald_rtol               = 1e-05
epsilon_r                = 6.7
ewald_geometry           = 3d

; OPTIONS FOR WEAK COUPLING ALGORITHMS = 
; Temperature coupling   = 
tcoupl                   = Berendsen
; Groups to couple separately = 
tc-grps                  = F Ca
; Time constant (ps) and reference temperature (K) = 
tau_t                    = 0.5 0.5
ref_t                    = 300.0 300.0
; Pressure coupling      = 
Pcoupl                   = Berendsen
Pcoupltype               = Isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar) = 
tau_p                    = 5.0
compressibility          = 1.20904E-06
ref_p                    = 1.0
; GENERATE VELOCITIES FOR STARTUP RUN = 
gen_vel                  = yes
gen_temp                 = 300.0
gen_seed                 = 173529
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

Any help is more than welcome.

Regards,
Adama

Adama Tandia
Modeling & Simulation
Corning INC
Corning, NY 14831 USA
Tel:  607 248 1036
Fax: 607 974 3405
www.corning.com