[CPMD-list] Possible bugs with NPT CPMD!
junwang
junwang at bigred.unl.edu
Fri Aug 4 23:42:53 CEST 2006
Dear All:
Recently I did a lot tests on finding a general procedure to implement NPT CPMD. But none of them was really successful.
First of all, I took one hydrogen molecule in a cubic box as the simplest test example(See the turtorial of CPMD day2.pdf).I followd the instructions to optimize wavefunctions first(same as inp1). Then I added several lines into inp2 to start a NPT CPMD simulation with temperature controlled by rescaling velocities.
INPUT:
&CPMD
MIRROR
MOLECULAR DYNAMICS
RESTART WAVEFUNCTION COORDINATES CELL
EMASS
200
TIMESTEP
2
TEMPCONTROL IONS
300 30
PARRINELLO-RAHMAN NPT
MAXSTEP
10000
TRAJECTORY SAMPLE
5000
&END
&SYSTEM
SYMMETRY
1
CELL
8.0 1.0 1.0 0.0 0.0 0.0
CUTOFF
30.
PRESSURE
0.0
ISOTROPIC CELL
&END
&DFT
FUNCTIONAL LDA
&END
&ATOMS
*H_MT_LDA.psp
LMAX=S
2
0.00000 0.00000 0.00000
0.80000 0.80000 0.80000
&END
OUTPUT:
****************************************************************
* AVERAGED QUANTITIES *
****************************************************************
MEAN VALUE <x> DEVIATION <x^2>-<x>^2
ELECTRON KINETIC ENERGY 0.455337E-04 0.148296E-04
IONIC TEMPERATURE 302.42 14.65
CELL TEMPERATURE 0.00 0.00
DENSITY FUNCTIONAL ENERGY -1.131362 0.477349E-03
CLASSICAL ENERGY -1.129926 0.466355E-03
CONSERVED ENERGY -1.129880 0.466823E-03
NOSE ENERGY ELECTRONS 0.000000 0.00000
NOSE ENERGY IONS 0.000000 0.00000
ION DISPLACEMENT 1.04871 0.669774
CELL VOLUME 512.000 0.00000
But the result showed that the cell size does not change at all during the simulation and the cell temperature is also zero all the time. I reported this problem to cpmd-list at cpmd.org. Juerg Hutter (hutter at pci.unizh.ch) told me that with a single molecule in a rather big box no forces on the box variables are induced, hence no movement of the box.
Then I began to use 16 atoms silicon system (the second example of day2.pdf tutorial) as my second test example. After the initial wavefunction optimization, I tried to start NPT simulation.
INPUT:
&CPMD
MIRROR
MOLECULAR DYNAMICS
RESTART WAVEFUNCTION COORDINATES CELL
EMASS
400
TIMESTEP
5
TEMPCONTROL IONS
1500 100
MAXSTEP
1000
PARRINELLO-RAHMAN NPT
SPLINE RANGE
5.0
&END
&DFT
FUNCTIONAL LDA
&END
&SYSTEM
SYMMETRY
FCC
SCALE S=2
CELL
20.526204 1.0 1.0 0 0 0
CUTOFF
8.0
PRESSURE
0.0
ISOTROPIC CELL
&END
&ATOMS
*Si_MT_LDA.psp KLEINMAN-BYLANDER
LMAX=P
16
0.00 0.00 0.00
0.25 0.25 0.25
1.00 0.00 0.00
1.25 0.25 0.25
0.00 1.00 0.00
0.25 1.25 0.25
1.00 1.00 0.00
1.25 1.25 0.25
0.00 0.00 1.00
0.25 0.25 1.25
1.00 0.00 1.00
1.25 0.25 1.25
0.00 1.00 1.00
0.25 1.25 1.25
1.00 1.00 1.00
1.25 1.25 1.25
&END
&BASIS
PSEUDO AO 2
0 1
&END
OUPUT:
NFI EKINC EKINH TEMPP EKS ECLASSIC EHAM DIS TCPU
1 8.76383 0.00000 1500.0 100.83262 100.93949 109.70332 0.652E-05 0.29
WARNING| SPLINE REGION SMALLER THAN MAXIMUM G-VALUE
2 47.52671 0.00000 1500.0 62.58080 62.68767 110.21438 0.322E-04 0.29
WARNING| SPLINE REGION SMALLER THAN MAXIMUM G-VALUE
3 77.97593 0.00000 1500.0 71.97072 72.07760 150.05353 0.865E-04 0.29
.. ...
Even if I increased SPLINE RANGE to 5.0, the warning" SPLINE REGION SMALLER THAN MAXIMUM G-VALUE" is still reported. In addition,some energy values are shown as "NaN".
If I controlled the temperature at 1500K by NOSE IONS, NOSE ELECTRONS first and start PARRINELLO-RAHMAN NPT later on, I still got the same error:
NFI EKINC EKINH TEMPP EKS ECLASSIC EHAM DIS TCPU
1 1.28022 NaN NaN -84.72569 NaN NaN 0.360E+10 0.09
WARNING| SPLINE REGION SMALLER THAN MAXIMUM G-VALUE
2 NaN NaN NaN NaN NaN NaN NaN 0.08
WARNING| SPLINE REGION SMALLER THAN MAXIMUM G-VALUE
3 NaN NaN NaN NaN NaN NaN NaN 0.09
I got a reply from Axel Kohlmeyer(akohlmey at cmm.chem.upenn.edu). He told me a few problems in my INPUT.
1. For variable cell you should provide a reference cell dimension that will be the largest possible cell for your dynamics.
2. Your plane wave cutoff is far too low. you have to make sure, that your stress tensor is converged wrt. to the basis set. So even if you get good forces (=structures with ~12 ry), you will need much more 40-50ry for the stress tensor. The way to decide if stress tensor is converged is to do single point calculations with increasing plane wave cutoff and print the stress tensor. After a while it should not change with cutoff.
He also told me that "NaN" may be related to the enforced symmetry 14. So he suggested me to redo the equilibration with symmetry 14.
Following his suggestions I increased the cutoff to 18, provided a reference cell and redid the equilibration with symmetry 14. Everything seemed normal during the equilibration. However, the error still occured when I turned on PARRINELLO-RAHMAN NPT.
INPUT:
** &CPMD **
** MIRROR **
** MOLECULAR DYNAMICS **
** RESTART WAVEFUNCTION COORDINATES **
** RESTART VELOCITIES NOSEP NOSEE **
** EMASS **
** 400 **
** TIMESTEP **
** 5 **
** NOSE IONS **
** 1500 1500 **
** NOSE ELECTRONS **
** 0.015 10000 **
** MAXSTEP **
** 2000 **
** TRAJECTORY SAMPLE **
** 200 **
** PARRINELLO-RAHMAN NPT **
** SPLINE RANGE **
** 5.0 **
** &END **
** **
** &DFT **
** FUNCTIONAL LDA **
** &END **
** &SYSTEM **
** SYMMETRY **
** 14 **
** SCALE S=2 **
** CELL **
** 20.526204 1.0 1.0 0 0 0 **
** REFERENCE CELL **
** 24.0 1.0 1.0 0 0 0 **
** CUTOFF **
** 18.0 **
** PRESSURE **
** 1.0 **
** ISOTROPIC CELL **
** &END **
OUTPUT:
NFI EKINC EKINH TEMPP EKS ECLASSIC EHAM DIS TCPU
1 1.29869 27.24130 393.5 -32.73910 -32.25976 -3.71978 0.153E+00 0.69
2 3.51333 NaN NaN -80.46287 NaN NaN 0.598E+07 0.70
WARNING| SPLINE REGION SMALLER THAN MAXIMUM G-VALUE
3 NaN NaN NaN NaN NaN NaN NaN 0.68
WARNING| SPLINE REGION SMALLER THAN MAXIMUM G-VALUE
4 NaN NaN NaN NaN NaN NaN NaN 0.69
Then he made a series of tests on one of his clusters and could confirm my problem. He found that there seems to be a bug with PARRINELLO-RAHMAN NPT in combination with the NOSE thermostat chains.If we use the regular PARRINELLO-RAHMAN it works fine.
Setting PARRINELLO-RAHMAN without NPT really solved the above problem, but the cell temperature is always zero and the cell length does not change.
INPUT:
** &CPMD **
** MIRROR **
** MOLECULAR DYNAMICS **
** RESTART WAVEFUNCTION COORDINATES **
** RESTART VELOCITIES NOSEP NOSEE **
** EMASS **
** 400 **
** TIMESTEP **
** 5 **
** NOSE IONS **
** 1500 1500 **
** NOSE ELECTRONS **
** 0.015 10000 **
** MAXSTEP **
** 1000 **
** TRAJECTORY SAMPLE **
** 100 **
** PARRINELLO-RAHMAN **
** SPLINE RANGE **
** 2.5 **
** STRESS TENSOR **
** 100 **
** &END **
OUTPUT:
NFI EKINC EKINH TEMPP EKS ECLASSIC EHAM DIS TCPU
1 0.02219 0.00000 1527.0 -62.43610 -62.32731 -62.30511 0.664E-05 0.65
2 0.01843 0.00000 1528.3 -62.43049 -62.32160 -62.30316 0.266E-04 0.66
3 0.01758 0.00000 1529.7 -62.42954 -62.32055 -62.30297 0.598E-04 0.64
4 0.02179 0.00000 1531.1 -62.43590 -62.32681 -62.30502 0.106E-03 0.63
5 0.01602 0.00000 1532.6 -62.42738 -62.31819 -62.30216 0.166E-03 0.65
.. ...
996 0.02079 0.00000 1865.0 -62.45735 -62.32447 -62.30368 0.242E+01 0.64
997 0.02080 0.00000 1864.9 -62.45735 -62.32448 -62.30368 0.243E+01 0.65
998 0.02078 0.00000 1864.8 -62.45732 -62.32446 -62.30368 0.243E+01 0.66
999 0.02070 0.00000 1864.6 -62.45723 -62.32438 -62.30368 0.243E+01 0.65
1000 0.02058 0.00000 1864.3 -62.45708 -62.32425 -62.30367 0.243E+01 0.65
.. ...
* AVERAGED QUANTITIES *
****************************************************************
MEAN VALUE <x> DEVIATION <x^2>-<x>^2
ELECTRON KINETIC ENERGY 0.194570E-01 0.756184E-03
IONIC TEMPERATURE 1616.18 135.59
CELL TEMPERATURE 0.00 0.00
DENSITY FUNCTIONAL ENERGY -62.438255 0.990662E-02
CLASSICAL ENERGY -62.323105 0.795040E-03
CONSERVED ENERGY -62.303648 0.100889E-03
NOSE ENERGY ELECTRONS 0.000000 0.00000
NOSE ENERGY IONS 0.000000 0.00000
ION DISPLACEMENT 1.10732 0.747528
CELL VOLUME 8647.81 0.391937
Then he replied to me that my cell volume is twice as large to what it should be for a bulk silicon crystal. In cases of a largely expanded cell, the 'force' on the cell are quite small. Finally he suggested me to start with a 64 atom primitive cubic cell like in the example in/out that was attached in his e-mail.
In his attached input file he did not specify "PRESSURE" or "ISOTROPIC CELL" in the $SYSTEM section.
INPUT:
&CPMD
MIRROR
MOLECULAR DYNAMICS CP
PARRINELLO RAHMAN
QUENCH BO ELECTRONS CELL
REAL SPACE WFN KEEP
TEMPCONTROL IONS
1500 200
EMASS
1000
TIMESTEP
4.0
SPLINE RANGE
5.0
TRAJECTORY SAMPLE
10
STRESS TENSOR
10
CONVERGENCE ORBITALS
1.0e-6
MAXSTEP
100
&END
&DFT
FUNCTIONAL LDA
&END
&SYSTEM
SYMMETRY
14
SCALE S=2
CELL
20.526204 1.0 1.0 0 0 0
REFERENCE CELL
22.00 1.0 1.0 0 0 0
CUTOFF
20.0
&END
&ATOMS
*Si_MT_LDA.psp KLEINMAN-BYLANDER
LMAX=P
64
0.50 1.00 0.50
1.00 0.50 0.50
0.50 0.50 0.00
1.00 1.00 0.00
0.75 0.25 0.75
0.25 0.25 0.25
0.25 0.75 0.75
0.75 0.75 0.25
1.50 1.00 0.50
2.00 0.50 0.50
1.50 0.50 0.00
2.00 1.00 0.00
1.75 0.25 0.75
1.25 0.25 0.25
1.25 0.75 0.75
1.75 0.75 0.25
0.50 2.00 0.50
1.00 1.50 0.50
0.50 1.50 0.00
1.00 2.00 0.00
0.75 1.25 0.75
0.25 1.25 0.25
0.25 1.75 0.75
0.75 1.75 0.25
1.50 2.00 0.50
2.00 1.50 0.50
1.50 1.50 0.00
2.00 2.00 0.00
1.75 1.25 0.75
1.25 1.25 0.25
1.25 1.75 0.75
1.75 1.75 0.25
0.50 1.00 1.50
1.00 0.50 1.50
0.50 0.50 1.00
1.00 1.00 1.00
0.75 0.25 1.75
0.25 0.25 1.25
0.25 0.75 1.75
0.75 0.75 1.25
1.50 1.00 1.50
2.00 0.50 1.50
1.50 0.50 1.00
2.00 1.00 1.00
1.75 0.25 1.75
1.25 0.25 1.25
1.25 0.75 1.75
1.75 0.75 1.25
0.50 2.00 1.50
1.00 1.50 1.50
0.50 1.50 1.00
1.00 2.00 1.00
0.75 1.25 1.75
0.25 1.25 1.25
0.25 1.75 1.75
0.75 1.75 1.25
1.50 2.00 1.50
2.00 1.50 1.50
1.50 1.50 1.00
2.00 2.00 1.00
1.75 1.25 1.75
1.25 1.25 1.25
1.25 1.75 1.75
1.75 1.75 1.25
&END
&BASIS
PSEUDO AO 2
0 1
&END
Using this INPUT I really got normal output as I expected.
OUTPUT:
NFI EKINC EKINH TEMPP EKS ECLASSIC EHAM DIS TCPU
1 0.00000 0.00038 1498.9 -252.65901 -252.21046 -252.21008 0.963E-05 1.97
2 0.00005 0.00146 1495.8 -252.65919 -252.21158 -252.21007 0.988E-04 2.00
3 0.00026 0.00302 1490.8 -252.65945 -252.21334 -252.21007 0.436E-03 2.01
4 0.00075 0.00479 1484.2 -252.65974 -252.21560 -252.21006 0.126E-02 2.02
5 0.00166 0.00648 1476.4 -252.65999 -252.21818 -252.21004 0.284E-02 2.01
.. ...
96 0.02158 0.00236 1321.2 -252.50792 -252.11257 -252.08862 0.165E+00 2.00
97 0.02221 0.00265 1309.9 -252.50546 -252.11348 -252.08862 0.176E+00 1.98
98 0.02268 0.00295 1500.0 -252.50283 -252.05396 -252.02834 0.187E+00 2.01
99 0.02300 0.00332 1487.6 -252.49981 -252.05466 -252.02833 0.199E+00 2.01
100 0.02321 0.00368 1475.1 -252.49663 -252.05522 -252.02833 0.212E+00 2.02
****************************************************************
* AVERAGED QUANTITIES *
****************************************************************
MEAN VALUE <x> DEVIATION <x^2>-<x>^2
ELECTRON KINETIC ENERGY 0.168993E-01 0.580536E-02
IONIC TEMPERATURE 1406.23 58.74
CELL TEMPERATURE 161.94 157.04
DENSITY FUNCTIONAL ENERGY -252.597747 0.495833E-01
CLASSICAL ENERGY -252.176940 0.514245E-01
CONSERVED ENERGY -252.157734 0.507740E-01
NOSE ENERGY ELECTRONS 0.000000 0.00000
NOSE ENERGY IONS 0.000000 0.00000
ION DISPLACEMENT 0.140194 0.753096E-01
CELL VOLUME 9016.36 168.223
When I specified the "PRESSURE" and "ISOTROPIC CELL", however, the kinetic energy of the cell became zero all the time and there is no change of cell size as shown below:
INPUT:
&SYSTEM
SYMMETRY
14
SCALE S=2
CELL
20.526204 1.0 1.0 0 0 0
REFERENCE CELL
22.00 1.0 1.0 0 0 0
CUTOFF
20.0
PRESSURE
0.0
ISOTROPIC CELL
&END
OUTPUT:
NFI EKINC EKINH TEMPP EKS ECLASSIC EHAM DIS TCPU
1 0.00000 0.00000 1499.2 -252.69024 -252.24161 -252.24160 0.438E-05 1.79
2 0.00005 0.00000 1497.1 -252.68963 -252.24164 -252.24159 0.175E-04 2.00
3 0.00024 0.00000 1493.6 -252.68877 -252.24182 -252.24158 0.394E-04 2.03
4 0.00068 0.00000 1489.0 -252.68783 -252.24225 -252.24157 0.700E-04 2.03
5 0.00143 0.00000 1483.7 -252.68698 -252.24300 -252.24157 0.109E-03 2.06
.. ...
96 0.00798 0.00000 1316.2 -252.52148 -252.12761 -252.11963 0.378E-01 2.08
97 0.00803 0.00000 1305.9 -252.51845 -252.12766 -252.11963 0.386E-01 1.79
98 0.00804 0.00000 1500.0 -252.51540 -252.06654 -252.05850 0.393E-01 1.80
99 0.00803 0.00000 1489.0 -252.51211 -252.06653 -252.05850 0.402E-01 1.80
100 0.00803 0.00000 1478.0 -252.50882 -252.06653 -252.05850 0.410E-01 2.02
****************************************************************
* AVERAGED QUANTITIES *
****************************************************************
MEAN VALUE <x> DEVIATION <x^2>-<x>^2
ELECTRON KINETIC ENERGY 0.787668E-02 0.222869E-02
IONIC TEMPERATURE 1410.45 57.92
CELL TEMPERATURE 0.00 0.00
DENSITY FUNCTIONAL ENERGY -252.624112 0.545877E-01
CLASSICAL ENERGY -252.202042 0.537170E-01
CONSERVED ENERGY -252.194166 0.535349E-01
NOSE ENERGY ELECTRONS 0.000000 0.00000
NOSE ENERGY IONS 0.000000 0.00000
ION DISPLACEMENT 0.139651E-01 0.123506E-01
CELL VOLUME 8648.22 0.125578E-01
CPU TIME 1.9812
He replied to me that there is a default value for pressure, 0kBar which is for all practical purposes.
Maybe the combination ISOTROPIC CELL/PARRINELLO-RAHMAN is is currently broken as well. He even doubted that I have changed the source code. But I never changed anything in the source codes.
Could anyone redo my tests and check if those problems are real bugs?
Best regards,
Sincerely,
Jun WANG
Department of Chemistry
University of Nebraska-Lincoln
536 Hamilton Hall
Lincoln, NE 68588-0304
U.S.A.
E-mail:junwang at unlserve.unl.edu
junwang at bigred.unl.edu
More information about the CPMD-list
mailing list