[CPMD-list] Energy components in PIMD

[Łukasz Walewski]

Dear All,

I have problems with total energy conservation in PIMD. All the EHAM
energy components are conserved except for ECLASSIC which systematically
drifts. The plots from my simulation are available at
http://bioexploratorium.pl/ljw/pimd/ I use CPMD v. 3.11.1 compiled with
IFC 10.1.008, ACML 4.0.1 and OpenMPI. Attached at the end of this e-mail
is the full input file.

According to my understanding the following holds:

ECLASSIC = EKS + EKIONS     (1)

On the other hand EKIONS should have something to do with the
temperature, however:

pi_md.F:420:              TEMPP(IPCURR)=EKINP*FACTEM*2.D0/GLIB_S
and
pi_md.F:418:              GLIB_S=3.D0*NAT

so this is EKINP that relates to the temperature. In my simulation EKINP
and so TEMPP are conserved. So this must be some component of ECLASSIC
that drifts. This is obviously not EKS (see the plot). So the suspect is
EKIONS, that I used to assume to be the kinetic energy of the ions. I
tried to track this down a little bit:

In the PIMD ECLASSIC is actually ECONSA (pi_md.F:583):

          IF(ENGPRI) WRITE(6,'(/,A,A)')
     *'    NFI  EKINC/P  TEMP   EKINP(PRI)  EKINP(VIR)      EKS/P ',
     *'     EQUANT    ECLASSIC        EHAM     TCPU'
          WRITE(6,'(I7,F8.5,F7.1,6F12.5,F9.2)')
     *      NFI,EKINA,TEMPA,QPKINP,QVKINP,ETOTA,EQUANT,ECONSA,EHAMA,TCPU
          WRITE(3,'(I7,F15.10,F7.1,6F20.10,F9.2)')
     *      NFI,EKINA,TEMPA,QPKINP,QVKINP,ETOTA,EQUANT,ECONSA,EHAMA,TCPU

Further on ECONSA becomes ECONS summed up over the PI replicas
(pi_md.F:576):

            ECONSA=ECONSA+ECONS(IP)

and finally ECONS itself becomes the sum of terms (pi_md.F:495):

            ECONS(IPCURR)=EKINP+ETOTV(IPCURR)/RNP+
     *                    ENOSE(IPCURR)/RNP+ENOSP(IPCURR)+
     *                    EHARV(IPCURR)+ECNSTR

from which I recognize only EKINP, ENOSE and ENOSP. I do not see how
this relates to (1).

Please note that I use termostat both for ions and electrons, so I would
expect kinetic energy to by conserved in such a job.

Could you please explain: 1) what is the meaning of the components of
ECLASSIC in PIMD, 2) what is the meaning of EKIONS ? and possibly give
some hints 3) why this value is not conserved in my simulation ?

With kind regards,
Lukasz


&CPMD
  PATH INTEGRAL
  MOLECULAR DYNAMICS CP
  RESTART WAVEFUNCTION COORDINATES LATEST
  TIMESTEP
    4.0
  EMASS
    400.0
  TEMPERATURE
    300.0
  NOSE ELECTRONS
    0.26 10000.0
  NOSE IONS MASSIVE
    300.0  180.8
  MAXSTEP
    2584
  SUBTRACT COMVEL ROTVEL
    10000
  TRAJECTORY XYZ SAMPLE
    1
  ISOLATED MOLECULE
  CENTER MOLECULE
  STORE
    1000
&END

&SYSTEM
  ANGSTROM
  SYMMETRY
    8
  CELL ABSOLUTE
    13.0 13.0 8.0 0.0 0.0 0.0
  CUTOFF
    30.0
&END

&DFT
  FUNCTIONAL BLYP
  GC-CUTOFF
    5.D-6
&END

&PIMD
  TROTTER DIMENSION
    16
  NORMAL MODES
    1.D0
  FACMASS
    1.D0
  DEBROGLIE CENTROID
    300.D0
&END

&ATOMS
*N_VDB_BLYP.psp FORMATTED
   LMAX=P
   4
       1.480434744000      -1.261526338000       0.012486679000
      -1.428836888000      -1.263268937000       0.020030520000
      -1.482423266000       1.260607946000       0.012029901000
       1.430161665000       1.261259358000       0.004480994000
*C_VDB_BLYP.psp FORMATTED
   LMAX=P
   20
       2.764321879000       1.617573818000      -0.000125061000
       3.855614684000       0.719804456000      -0.000057419000
       3.883340408000      -0.678014026000       0.004378686000
       2.816852951000      -1.603192445000       0.010125751000
       2.919542448000      -3.046779447000       0.014515219000
       1.637225743000      -3.547145051000       0.019451505000
       0.739381956000      -2.412298909000       0.018118048000
      -0.677939805000      -2.411241253000       0.021787063000
      -1.565841060000      -3.552172948000       0.027767194000
      -2.851571375000      -3.063067365000       0.029521295000
      -2.763080194000      -1.618905630000       0.024634604000
      -3.853925457000      -0.720052370000       0.024562307000
      -3.883570717000       0.678208069000       0.020129751000
      -2.818027537000       1.604764172000       0.014380689000
      -2.917826697000       3.049092800000       0.009981340000
      -1.634343765000       3.546387751000       0.005051935000
      -0.738698022000       2.409860388000       0.006396510000
       0.678687181000       2.408624576000       0.002727905000
       1.566148163000       3.550017372000      -0.003252572000
       2.852089024000       3.061544408000      -0.005009259000
*H_VDB_BLYP.psp FORMATTED
   LMAX=S
   14
       4.829896620000       1.209449328000      -0.004160981000
       4.874865030000      -1.131837359000       0.003273118000
       3.847285348000      -3.609138791000       0.013924943000
       1.339951403000      -4.589537686000       0.023583475000
      -1.257219301000      -4.591087509000       0.030318008000
      -3.775889801000      -3.630740840000       0.033747037000
      -4.828126595000      -1.209793447000       0.028665970000
      -4.875902354000       1.130204086000       0.021243303000
      -3.843471723000       3.614955577000       0.010554882000
      -1.334243713000       4.588096584000       0.000914849000
       1.257025926000       4.588704910000      -0.005801358000
       3.776400021000       3.629099090000      -0.009234599000
       1.216317908000      -0.004972921000       0.009118666000
      -1.218997172000       0.003787130000       0.015400566000
&END