[gmx-users] Problems with implementation of tabulated potentials of Buckingham function

[Diego Andres Henao Gonzalez]

Hello everyone.

I'm currently simulating monoatomic argon at a pressure of 20 bar and a
temperature of 83 kelvin. For information reported in the literature, under
these conditions argon behaves as a liquid and its density is approximately
1400 kg / m³. This system works with a Buckingham potential. I have run the
simulations (minimization, NVT at 200 ps and NPT at 1000 ps) to calculate
the density using the internal code of gromacs for this potential by
specifying the following topology:

[ defaults ]

; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ

2 1 no 1.0 1.0

#include "./atomtypes.itp"

#include "./ar.itp"


[ system ]

LIG

[ molecules ]

LIG 343


The files atomtypes.itp and ar.itp are the following:


*for atomtypes.itp:*

[ atomtypes ]

;name at.num mass charge ptype a(kJ/mol) b(nm^-1) c6(kJ/mol*nm^6)

; Argon gas

AR 18 39.950 0 A 1017718.0 36.6 0.006142


*for ar.itp:*

[ moleculetype ]

; molname nrexcl

LIG 1

[ atoms ]

; id attype resnr residuname atname cgnr charge mass atB qB mB

1 AR 1 AR AR 1 0.0 39.95000 AR 0.0 39.95000 ;


In these conditions, the simulation works well and does not generate
errors, giving a density according to the reported information.

The problem that I have is related to the use of tables for Buckingham's
potential. I have tried many types of combinations for start the
simulation, but in all cases I get the following error:

*NOTE: The number of threads is not equal to the number of (logical) cores*

*and the -pin option is set to auto: will not pin thread to cores.*

*This can lead to significant performance degradation.*

*Consider using -pin on (and -pinoffset in case you run multiple jobs).*


*starting mdrun 'LIG'*

*100000 steps, 200.0 ps.*

*Segment violation (`core 'generated)*



I have built the tables in terms of Buckingham's potential as follows:

1) The exponential function only depending on the radius: Exp (-r) and
attractive term in terms of radius
2) The exponential function depending on the radius and the parameter B:
Exp (-r * B) and attractive term in terms of the radius
3) The exponential function as a function of the radius, the parameter B
and the parameter A: A * Exp (-r * B) and attractive term in terms of the
radius
4) The exponential function as a function of radius, parameter B and
parameter A: A * Exp (-r * B) and attractive term in terms of radius and
parameter C


I have configured each of these combinations in different ways in the
[default] and I have placed the parameters A B C them as appropriate (for
example, if the table is for case 2, I don't define the parameter B in the
topology).

For all cases I used the information reported in the following link:
http://www.gromacs.org/Documentation/How-tos/Tabulated_Potentials

However, after many attempts, the simulation generates the same error. I
have the following questions:

1) How should I build the tables for Buckingham? (ejm: Exp (-r) or Exp (-r
* B) or A * Exp (-r * B))
2) How should I define parameters A B and C of Buckingham in the
[atomtypes] and in what order should I place them? (ejm: I define only A
and C, or I put A B and C)
3) How is the combination rule defined in the [defaults]? I put 1 or 2 and
if I put some as I define the atomtypes?


Finally, since there is only one interaction, I'm defining a single table
in the .mdp file:

vdwtype = User
energygrps = AR

Is the foregoing well defined?

I'm going to leave all the files to work this simulation:


*ar.itp:*


; Topology file for ARGON model

[ moleculetype ]

; molname nrexcl

LIG 1

[ atoms ]

; id attype resnr residuname atname cgnr charge mass atB qB mB gauss

1 AR 1 AR AR 1 0.0 39.95000 AR 0.0 39.95000



*atomtypes.itp*

[ atomtypes ]

;name at.num mass charge ptype a(kJ/mol) b(nm^-1) c6(kJ/mol*nm^6)

; Argon Noble gas

AR 18 39.950 0 A 1.0 1.0 1.0


*argon.gro*

1

1LIG AR 1 0.231 -0.387 -0.000

0.00000 0.00000 0.00000


*table.cpp*


#include <stdio.h>

#include <math.h>


main()

{

FILE *fout;

double r;


/* double B;*/

/* B=35.6;*/

double A_A;

double C_A;

double deltar;

double B_A;

A_A = 1017718.0;

B_A = 36.6;

C_A = 0.00614244;

deltar = 0.000001;

fout = fopen("table_AR_AR.xvg", "w");

fprintf(fout, "#\n# Buckingham Potential_ONLY VDW\n#\n");


for (r=0; r<=2; r+=0.0005) {


/* f es electrostaticas; g y h vdw*/


double f = 0;

double fprime = 0;

double g = ((-1)/(pow(r,6));

double gmas = (-1)/(pow((r+deltar),6));

double gmen = ((-1))/(pow((r-deltar),6));

double gprime = (-1)*(gmas - gmen)/(2*deltar);


double h = exp((-r));


double hmas = exp((-r)+deltar);

double hmen = exp((-r)-deltar);

double hprime = (-1)*(hmas - hmen)/(2*deltar);


/* print output */

if (r<0.04) {

fprintf(fout, "%12.10e %12.10e %12.10e %12.10e %12.10e %12.10e %12.10e\n",
r,0.0,0.0,0.0,0.0,0.0,0.0);

} else {

fprintf(fout, "%12.10e %12.10e %12.10e %12.10e %12.10e %12.10e %12.10e\n",
r,f,fprime,g,gprime,h,hprime);

}

}


fclose(fout);

return(0);

}


*topol.top*

[ defaults ]

; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ

2 1 no 1.0 1.0

#include "./atomtypes.itp"

#include "./ar.itp"

[ system ]

LIG

[ molecules ]

LIG 343


*nvt.mdp*

; Run control

integrator = sd ; Langevin dynamics

tinit = 0

dt = 0.002

nsteps = 100000 ; 100 ps

nstcomm = 100

; Output control

nstxout = 500

nstvout = 500

nstfout = 0

nstlog = 500

nstenergy = 500

nstxout-compressed = 0

; Neighborsearching and short-range nonbonded interactions

cutoff-scheme = group

nstlist = 20

ns_type = grid

pbc = xyz

rlist = 1.0

; Electrostatics

coulombtype = PME

rcoulomb = 1.0

; van der Waals

vdwtype = User

vdw-modifier = none

rvdw-switch = 0.9

rvdw = 1.0

; Apply long range dispersion corrections for Energy and Pressure

DispCorr = EnerPres


; Extension of the potential lookup tables beyond the cut-off

table-extension = 1

; Separate tables between energy group pairs

energygrps = AR


; Spacing for the PME/PPPM FFT grid

fourierspacing = 0.12

; EWALD/PME/PPPM parameters

pme_order = 6

ewald_rtol = 1e-06

epsilon_surface = 0

; Temperature coupling

; tcoupl is implicitly handled by the sd integrator

tc_grps = system

tau_t = 1.0

ref_t = 300

; Pressure coupling is off for NVT

Pcoupl = no

tau_p = 0.5

compressibility = 4.5e-05

ref_p = 1.0

; Generate velocities to start

gen_vel = yes

gen_temp = 300

gen_seed = -1

; options for bonds

constraints = all-bonds

; Type of constraint algorithm

constraint-algorithm = lincs

; Do not constrain the starting configuration

continuation = no

; Highest order in the expansion of the constraint coupling matrix

lincs-order = 12


Thank you very much. Regards



-- 
I.Q. Diego Andres Henao Gonzalez
M.sc Candidate
Department of Chemical and Environmental Engineering (Bogotá)
National University of Colombia