[CPMD-list] Bug in center of mass velocity calculation with fixed atoms

[Ata Roudgar]

Hi everybody,

I just found a bug in CPMD and I thought I have to report it.
I am using CPMD version 3.11.1.


The bug is related to the calculation of center of mass velocity in the case 
when the constraints are used in ATOM section. So if you do not use any 
constraint this bug does not affect on you. 

The center of mass velocity is calculated in subroutine RINVEL at the 
beginning of each MD where the initial velocities are calculated according to 
maxwell Boltzmann distribution or random shooting. The center of mass 
velocity calculation is necessary in order to avoid the motion of the center 
of the mass of the system during the MD simulation.

The bug appears when constrains are in used. In this case the center of mass 
is calculated according to the following formula


V_{cm}(i) = \sigma_{j,k}{Vel(i,j,k)*m_{j,k}} / M

where i is (x,y,z) and j and j and k are atom species and number of atom per 
species respectively. M is the total mass. Later on the Vel is define as 
following:

Vel =Vel -Vcm

in this manner the center of the mass of the system will be always zero. Note 
that those Vel which are listed in the constains will be consider as zero in 
SUBROUTINE   TAUCL.  The total mass M is given by:

M = sigma_{i,j}{m_{i,j}}
where i and j are index for atom species and number of atom per species BUT 
m_{i,j} MUST be zero for those atoms which are listed in the CONSTRAINS 
section and are supposed to be fixed. Because when we do not consider any 
velocity for fixed atoms we do not consider any mass for them too otherwise 
we always receive the following message:

 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !! RINVEL| WARNING! CENTER OF MASS VELOCITY NOT ZERO          !!
 !!    0.0000024859   0.0000006602   0.0000026434              !!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



Here the correction:
File rinvel.F


line 15:
      INCLUDE 'cotr.inc'
line 27 add:
      INTEGER    IAA
line 113 add:
      IAA=0
line 116 add:
      IAA=IAA+1
line 120 replace:
PMA00=PMA00+PMA0(IS) -----> 

PMA00=PMA00+(1.0/3.0)*(LSKCOR(1,IAA)+
     &               LSKCOR(2,IAA)+LSKCOR(3,IAA))*PMA0(IS)

      

line 120:
              PMA00=PMA00+PMA0(IS)


 line 161 add:
      IAA=0
line 164 add:
      IAA=IAA+1
line 167 replace:
PMA00=PMA00+PMA0(IS) -----> 

PMA00=PMA00+(1.0/3.0)*(LSKCOR(1,IAA)+
     &               LSKCOR(2,IAA)+LSKCOR(3,IAA))*PMA0(IS)




I tested the program and now it works fine,
Thanks for your time,
Ata