[gmx-users] Problems with implementation of tabulated potentials of Buckingham function
Diego Andres Henao Gonzalez
dhenaog at unal.edu.co
Wed Apr 25 22:04:49 CEST 2018
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
More information about the gromacs.org_gmx-users
mailing list