[gmx-users] Strange behaviour of added atomtype in simulation

[Kamps, M.]

Dear gmx-user.

I am trying to implement an an-organic ionic-bound surface into GROMACS.
I've added the residue (which is Al, aluminium, uncharged) into an rtp
file, and added this to an edited Charmm force-field. The atom is added to
atomtypes.atp and ffnonbonded.itp. The LJ potentials are found from the
literature.
The surface is created by insert-molecules the individual atoms in a
specific FCC orientation. I then simulated the surface in a vacuum to see
how it behaves. To enable this, the surface is energy minimized and NVT and
NPT equilibriated.

The problem arises when it is simulated at room temperature (300K), where
the surface will behave like a liquid, instead of a solid. This seems
strange since the critical temperature is not reached by far. I tried
testing the surface at a lower temperature (10K) and there it behaves as
expected. At 10K the surface will stay in shape, in its optimal FCC
configuration and it will slightly vibrate as expected.
If the temperature is set to 300K, the surface will immediately deform and
atoms will move all through each other and around each other. The surface
will keep its shape, due to atomistic interactions, however it will not
stay in its optimal FCC configuration, which is what should happen.

What am I doing wrong?

The used md.mdp file is given below:
title = Al Surface MD simulation
; Run parameters
integrator = md-vv ; leap-frog integrator
nsteps = 200000 ; 1 * 200000 = 200 ps (0.2 ns)
dt     = 0.001 ; 1 fs
; Output control
nstxout       = 1000 ; save coordinates every 1.0 ps
nstvout         = 1000 ; save velocities every 1.0 ps
nstenergy         = 1000 ; save energies every 1.0 ps
nstlog         = 1000 ; update log file every 1.0 ps
; Bond parameters
continuation       = yes ; Do not continue
constraint_algorithm     = lincs     ; holonomic constraints
constraints             = all-bonds ; all bonds (even heavy atom-H bonds)
constrained
lincs_iter             = 1     ; accuracy of LINCS
lincs_order             = 4     ; also related to accuracy
; Neighborsearching
cutoff-scheme   = Verlet
ns_type     = grid ; search neighboring grid cells
nstlist     = 10     ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb     = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw     = 1.0 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype     = cut-off ; Particle Mesh Ewald for long-range
electrostatics
; Temperature coupling is on
tcoupl = berendsen        ; modified Berendsen thermostat
tc-grps = AL ; one coupling group
tau_t = 0.1          ; time constant, in ps
ref_t = 300          ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl       = berendsen ; Pressure coupling on in NPT
pcoupltype         = semiisotropic        ; scaling with different value in
z direction
tau_p         = 2.0        ; time constant, in ps
ref_p         = 0.0 0.0        ; reference pressure, for both directions,
in bar
compressibility     = 0.0 4.5e-5        ; compressibility, for both
directions, in bar^-1 (4.5e-5)
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = no ; velocity generation is off
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme