Help please: constrained MD does not work

[5-903]

Hello. I would like to perform regular constrained MD, fixing a selected interatomic distance. Following the instructions in the manual, I've created a file named ICONST and specified the distance to be fixed, i.e.:

R 4 20 0

(note this is the **entire** content of the file -- no comment or blank lines)

I've run MD with that ICONST file in my working folder (see the INCAR and POSCAR files below), but the above specified constraint appears not to have any effect at all. Visualization (vasprun.xml file by VMD) reveals that the selected distance (an OH bond) shrinks and elongates, as it does in a regular MD without constraints. No REPORT file was produced (although it should be there, according to the manual).

I've tried this with two versions of VASP (5.2.12 and 5.3.5) and the result is the same in both cases -- the program appears to be completely ignoring the ICONST file, running regular unconstrained MD. Also, I found no mention of constraints in the OUTCAR file, or any indication that the ICONST file has been considered in the run.

What did I possibly do wrong? Is there something to be added to the ICONST or INCAR file?

I would appreciate any comments and suggestions -- thanks in advance!

*** INCAR file: ***
LWAVE = .TRUE. do not write WAVECAR file
LCHARG = .TRUE. do not write CHGCAR file
LVTOT = .FALSE. do not write LOCPOT file
LELF = .FALSE.

NPAR = 2

Electronic Relaxation 1
PREC = Medium Medium=default, Low, High; affects ENMAX, mesh, pspot
NELM = 100 max number of electronic steps
EDIFF = 1E-04 energy stopping-criterion for electr. iterations
EDIFFG = -.003 force stopping-criterion for geometry steps
NSIM = 40 stores dielectric matrix longer and improves convergence

Ionic Relaxation
NSW = 30000 max number of geometry steps
IBRION = 0 ionic relax: 0-MD 1-quasi-New 2-CG 3-damped dynamics
ISIF = 0 (0:force=y stress=n ions=y shape=n volume=n,
ISYM = 0 1=use symmetry, 0 = no symmetry
POTIM = 1.0 time step for geo-opt (reduce if small distances exist) and MD (in fs)

MDALGO = 2
SMASS = 0.5
TEBEG = 300.0

DOS related values:
ISMEAR = 1 -4-tet -1-fermi 1=Methfessel/Paxton 1.order
SIGMA = .05 broadening in eV

Electronic Relaxation
IALGO = 48 algorithm (8=CG for small, 48=RMM for big systems)
LREAL = A real-space projection

*** POSCAR file: ***
Generated by cif2cell 1.1.5 from CSD reference: OXACDH26. C2 H6 O6 : D.Zobel et al., Acta Crystallogr.,Sect.B:Struct.Sci. 48, 837- (1992). Species order: H C O
6.093000
1.000000000000000 0.000000000000000 0.000000000000000
0.000000000000000 0.569341867717052 0.000000000000000
-0.529324970666470 0.000000000000000 1.884396002353298
12 4 12
Direct
0.029000000000000 0.008000000000000 0.214000000000000
0.971000000000000 0.992000000000000 0.786000000000000
0.471000000000000 0.508000000000000 0.286000000000000
0.529000000000000 0.492000000000000 0.714000000000000
0.442000000000000 0.691000000000000 0.119000000000000
0.558000000000000 0.309000000000000 0.881000000000000
0.058000000000000 0.191000000000000 0.381000000000000
0.942000000000000 0.809000000000000 0.619000000000000
0.636000000000000 0.470000000000000 0.150000000000000
0.364000000000000 0.530000000000000 0.850000000000000
0.864000000000000 0.970000000000000 0.350000000000000
0.136000000000000 0.030000000000000 0.650000000000000
0.955170000000000 0.059520000000000 0.052080000000000
0.044830000000000 0.940480000000000 0.947920000000000
0.544830000000000 0.559520000000000 0.447920000000000
0.455170000000000 0.440480000000000 0.552080000000000
0.085280000000000 0.944620000000000 0.150360000000000
0.914720000000000 0.055380000000000 0.849640000000000
0.414720000000000 0.444620000000000 0.349640000000000
0.585280000000000 0.555380000000000 0.650360000000000
0.778510000000000 0.244980000000000 0.036320000000000
0.221490000000000 0.755020000000000 0.963680000000000
0.721490000000000 0.744980000000000 0.463680000000000
0.278510000000000 0.255020000000000 0.536320000000000
0.548390000000000 0.634770000000000 0.178460000000000
0.451610000000000 0.365230000000000 0.821540000000000
0.951610000000000 0.134770000000000 0.321540000000000
0.048390000000000 0.865230000000000 0.678460000000000

[admin]

To perform constrained MD "vasptbdyn" executable should be started. Constrains are indicated in REPORT as ">Const_coord"

[5-903]

Thank you for pointing me to the tbdyn module. Apparently I have not included it at the preprocessing stage. In fact I have invoked the -Dtbdyn flag in the CPP options, recompiled the code, and now constrained dynamics seems to be working normally with the vasp binary (at least according to preliminary tests). The slow growth approach also works.

I have another question though. How can I monitor the force acting along the constrained variable? I need the force and its dynamics in order to perform thermodynamic integration. In the REPORT file there is some data about energetics (see below). Is any of these quantities (E_const?) related to the dynamical force acting along the constraint or should I derive the force solely from data in the OUTCAR file? What is the exact meaning of "E_const", "EPS" and "ES" anyway? Unfortunately, the content of the REPORT file appears to be poorly documented in the manual, stating only that it includes quantities needed to compute the free-energy gradient (which is probably what I need), but lacking detailed description.

>Energies
E_tot E_pot E_kin E_const EPS ES
e_b> -0.15597290E+03 -0.15992493E+03 0.39390922E+01 -0.19042448E-02 0.10611393E-01 0.42287172E-02

Many thanks for your help again!

Jernej

[admin]

The force you are interested in (=gradient of the free energy) is in REPORT printed as "lambda" (the number behind b_m>).
Have a look also at corresponding equations:
http://cms.mpi.univie.ac.at/vasp/vasp/C ... amics.html

[5-903]

Hmmm... no trace of "lambda" or "b_m>" in my REPORT file (see below -- a full listing of a step is given). Any idea how to invoke the printing?

Jernej

========================================
MD step No. 2
========================================

SHAKE converged in 21 steps

>Const_coord
cc> R 1.08954 1.08954 -0.339476E-08

>Energies
E_tot E_pot E_kin E_const EPS ES
e_b> -0.15584101E+03 -0.16296839E+03 0.66829433E+00 -0.11925842E-04 0.66336500E-01 0.63927585E+01

>Temperature
T_sim T_inst
tmprt> 300.000 193.880

>Thermostat:
s(n+1) sdot(n+1) s(n) sdot(n) Qs
s_b> 3.08252 -0.00870 3.09102 -0.00830 0.50000

RANDOM_SEED = 411848287 411848324

[admin]

The output of the blue-moon sampling is switched on with LBLUEOUT=.TRUE.

[5-903]

Thank you, with LBLUEOUT this works as you suggested.

I'd appreciate further clarification though. In my case it looks like the application of eq. (6.29) (mathematical derivation of free energy gradients), yields very simple form, namely (dA/dksi) = <lambda>. This is due to the fact that the |Z|^(-1/2) remains constant throughout the constrained dynamics (fixed interatiomc distance), and that the lenghty term added to lambda in the expression (6.29) ("GkT"?) seems to be negligible (a typical example is given below). The entire expression then simplifies to the value of lambda.

>Blue_moon
lambda |z|^(-1/2) GkT |z|^(-1/2)*(lambda+GkT)
b_m> -0.181209E+02 0.970143E+00 -0.450497E-14 -0.175799E+02

HOWEVER, the resulting values of lambda appear not to be consistent with the "dumb" test, made by computing the net force acting on the involved two atoms, projected onto the constrained vector. This force is derivable from the list of forces in the OUTCAR file. It looks like the net force on constraint, as computed from the OUTCAR file, is approximately twice as large as the value of lambda. However, much worse is the fact that the dynamical profile of this force is not identical to that of lambda.

Regardless of my double-checking, I wonder what the units of lambda are in the case of constrained distance. I would assume lambda should be given in eV/(angstrom)^2. Correct?

[admin]

The unit of lambda depends on the constrained coordinate. When constrained coordinate is Angstroem,
then the unit is eV/A. The integration of the gradient of the free energy (last number in the b_M> line)
along Angstroems (eq. 6.30) then gives the change of the free energy in eV.