TABLE OF CONTENTS


ABINIT/m_pimd [ Modules ]

[ Top ] [ Modules ]

NAME

  m_pimd

FUNCTION

  This module provides several routines and datatypes for the
  Path-Integral Molecular Dynamics (PIMD) implementation.

COPYRIGHT

 Copyright (C) 2010-2022 ABINIT group (GG,MT)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 MODULE m_pimd
24 
25  use defs_basis
26  use m_abicore
27  use m_dtset
28  use m_errors
29  use m_io_tools
30  use m_random_zbq
31 
32  use m_numeric_tools,  only : uniformrandom
33  use m_symtk,          only : matr3inv
34  use m_geometry,       only : mkradim
35 
36  implicit none
37 
38  private
39 
40 !public procedures
41  public :: pimd_init
42  public :: pimd_nullify
43  public :: pimd_destroy
44  public :: pimd_init_qtb
45  public :: pimd_skip_qtb
46  public :: pimd_print
47  public :: pimd_is_restart
48  public :: pimd_temperature
49  public :: pimd_initvel
50  public :: pimd_langevin_random
51  public :: pimd_langevin_random_qtb
52  public :: pimd_langevin_random_bar
53  public :: pimd_langevin_random_init
54  public :: pimd_energies
55  public :: pimd_forces
56  public :: pimd_langevin_forces
57  public :: pimd_nosehoover_forces
58  public :: pimd_stresses
59  public :: pimd_diff_stress
60  public :: pimd_predict_taylor
61  public :: pimd_predict_verlet
62  public :: pimd_nosehoover_propagate
63  public :: pimd_coord_transform
64  public :: pimd_force_transform
65  public :: pimd_apply_constraint
66  public :: pimd_mass_spring

m_pimd/pimd_apply_constraint [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_apply_constraint

FUNCTION

  Modify forces to take into account an holonomic constraint
  according to "pimd_constraint" parameter
  Available constraints:
    0: no constraint
    1: linear combination of coordinates

INPUTS

  constraint=type of constraint to be applied
  mass(natom)=fictitious masses of atoms
  natom=number of atoms
  trotter=Trotter number
  wtatcon(3,natom)=weights for atomic constraints
  xcart(3,natom,trotter)=cartesian coordinates of atoms

OUTPUT

  constraint_output(2)=several (real) data to be output
                       when a constraint has been applied

SIDE EFFECTS

  forces(3,natom,trotter)=array containing forces

SOURCE

2367 subroutine pimd_apply_constraint(constraint,constraint_output,forces,mass,natom,&
2368 &                                trotter,wtatcon,xcart)
2369 
2370 !Arguments ------------------------------------
2371 !scalars
2372  integer,intent(in) :: constraint,natom,trotter
2373 !arrays
2374  real(dp),intent(in) :: mass(natom),wtatcon(3,natom),xcart(3,natom,trotter)
2375  real(dp),intent(out) :: constraint_output(2)
2376  real(dp),intent(inout) :: forces(3,natom,trotter)
2377 !Local variables-------------------------------
2378 !scalars
2379  integer :: iatom,ii,iimage
2380  real(dp) :: af,lambda_cst,masstot,one_over_trotter,xcart_centroid,zz
2381  !character(len=500) :: msg
2382 !arrays
2383  real(dp) :: force_centroid(3),lambda_com(3),mat(3,3),matinv(3,3),vec(3),weightsum(3)
2384 
2385 !************************************************************************
2386 
2387 
2388  select case(constraint)
2389 
2390 !=== No constraint =======================================================
2391  case(0)
2392 
2393    constraint_output(:)=zero
2394    return
2395 
2396 !=== Linear combination of centroid coordinates ==========================
2397  case(1)
2398 
2399 !  Some useful quantities
2400    masstot=sum(mass)*dble(trotter)
2401    weightsum(:)=sum(wtatcon,dim=2)*dble(trotter)
2402    zz=zero;af=zero;force_centroid=zero
2403    do iimage=1,trotter
2404      do iatom=1,natom
2405        do ii=1,3
2406          force_centroid(ii)=force_centroid(ii)+forces(ii,iatom,iimage)
2407          af=af+wtatcon(ii,iatom)*forces(ii,iatom,iimage)/mass(iatom)
2408          zz=zz+wtatcon(ii,iatom)**2/mass(iatom)
2409        end do
2410      end do
2411    end do
2412    vec(:)=force_centroid(:)-weightsum(:)*af/zz
2413    do ii=1,3
2414      mat(:,ii)=-weightsum(:)*weightsum(ii)/zz
2415      mat(ii,ii)=mat(ii,ii)+masstot
2416    end do
2417    call matr3inv(mat,matinv)
2418 
2419    !Calculation of a Lagrange multipliers:
2420    ! lambda_cst: to apply the constraint
2421    ! lambda_com: to maintain the position of the center of mass
2422    lambda_com(:)=matmul(matinv,vec)
2423    lambda_cst=(af-dot_product(weightsum,lambda_com))*dble(trotter)/zz
2424 
2425    !Modification of forces
2426    one_over_trotter=one/dble(trotter)
2427    do iimage=1,trotter
2428      do iatom=1,natom
2429        do ii=1,3
2430          forces(ii,iatom,iimage)=forces(ii,iatom,iimage) &
2431 &                               -lambda_cst*wtatcon(ii,iatom)*one_over_trotter &
2432 &                               -lambda_com(ii)*mass(iatom)
2433        end do
2434      end do
2435    end do
2436 
2437    !Computation of relevant outputs
2438    constraint_output(:)=zero
2439    !1-Reaction coordinate
2440    do iatom=1,natom
2441      do ii=1,3
2442        xcart_centroid=zero
2443        do iimage=1,trotter
2444          xcart_centroid=xcart_centroid+xcart(ii,iatom,iimage)
2445        end do
2446        constraint_output(1)=constraint_output(1)+xcart_centroid*wtatcon(ii,iatom)
2447      end do
2448    end do
2449    constraint_output(1)=constraint_output(1)/dble(trotter)
2450    !2-Force on reaction coordinate
2451    constraint_output(2)=-lambda_cst
2452 
2453  end select
2454 
2455 end subroutine pimd_apply_constraint

m_pimd/pimd_coord_transform [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_coord_transform

FUNCTION

  Apply a coordinate transformation on a given vector field
  (defined for each atom in in cell) - for PIMD
  Possible choices for the transformation:
    0: no transformation
    1: normal mode tranformation
    2: staging transformation

INPUTS

  ioption=option given the direction of the transformation
     +1: from primitive coordinates to transformed coordinates
     -1: from transformed coordinates to primitive coordinates
  natom=number of atoms
  transform=coordinate transformation:
            0: no tranformation
            1: normal mode transformation
            2: staging transformation
  trotter=Trotter number

OUTPUT

SIDE EFFECTS

  array(3,natom,trotter)=array to be transformed

SOURCE

2064 subroutine pimd_coord_transform(array,ioption,natom,transform,trotter)
2065 
2066 !Arguments ------------------------------------
2067 !scalars
2068  integer,intent(in) :: ioption,natom,transform,trotter
2069 !arrays
2070  real(dp),intent(inout) :: array(3,natom,trotter)
2071 !Local variables-------------------------------
2072 !scalars
2073  integer :: iatom,ii,iimage,iimagem,iimagep,iimage2,ll
2074 !arrays
2075  real(dp),allocatable :: array_temp(:,:,:)
2076  real(dp),allocatable :: nrm(:,:)
2077 
2078 !************************************************************************
2079 
2080 !=== No transformation ===================================================
2081  if (transform==0) then
2082 
2083    return
2084 
2085 !=== Normal mode transformation ==========================================
2086  else if (transform==1) then
2087 
2088 !  ---------------- From primitive to transformed coordinates ------------
2089    if (ioption==+1) then
2090      ABI_MALLOC(array_temp,(3,natom,trotter))
2091 
2092      array_temp(:,:,1)=zero
2093      do iimage=1,trotter
2094        array_temp(:,:,1)=array_temp(:,:,1)+array(:,:,iimage)
2095      end do
2096 
2097      do ll=2,trotter/2
2098        array_temp(:,:,2*ll-2)=zero
2099        array_temp(:,:,2*ll-1)=zero
2100        do iimage=1,trotter
2101          array_temp(:,:,2*ll-2)=array_temp(:,:,2*ll-2)+array(:,:,iimage)* &
2102 &                              cos(two_pi*dble((iimage-1)*(ll-1))/dble(trotter))
2103          array_temp(:,:,2*ll-1)=array_temp(:,:,2*ll-1)-array(:,:,iimage)* &
2104 &                              sin(two_pi*dble((iimage-1)*(ll-1))/dble(trotter))
2105        end do
2106      end do
2107 
2108      array_temp(:,:,trotter)=zero
2109      do iimage=1,trotter
2110        array_temp(:,:,trotter)=array_temp(:,:,trotter)+&
2111 &                              array(:,:,iimage)*dble((-1)**(iimage-1))
2112      end do
2113 
2114      array=array_temp/dble(trotter)
2115 
2116      ABI_FREE(array_temp)
2117 
2118 !  ---------------- From transformed to primitive coordinates ------------
2119    else if (ioption==-1) then
2120 
2121     ABI_MALLOC(array_temp,(3,natom,trotter))  !real part
2122     ABI_MALLOC(nrm,(trotter,trotter))  !real part
2123 
2124    do iimage=1,trotter
2125      nrm(iimage,1)=one
2126      nrm(iimage,trotter)=dble((-1)**(iimage-1))
2127    end do
2128 
2129    do ll=2,trotter/2
2130      do iimage=1,trotter
2131        nrm(iimage,2*ll-2)= two*cos(two_pi*dble((ll-1)* &
2132 &                        (iimage-1))/dble(trotter))
2133        nrm(iimage,2*ll-1)=-two*sin(two_pi*dble((ll-1)* &
2134 &                        (iimage-1))/dble(trotter))
2135      end do
2136    end do
2137 
2138     do iimage=1,trotter
2139       array_temp(:,:,iimage)=zero
2140       do ll=1,trotter
2141         array_temp(:,:,iimage)=array_temp(:,:,iimage)+nrm(iimage,ll)*array(:,:,ll)
2142       end do
2143     end do
2144     array=array_temp
2145 
2146     ABI_FREE(array_temp)
2147     ABI_FREE(nrm)
2148 
2149    end if ! ioption
2150 
2151 !=== Staging transformation ==============================================
2152  else if (transform==2) then
2153 
2154 !  ---------------- From primitive to transformed coordinates ------------
2155    if (ioption==+1) then
2156 
2157      ABI_MALLOC(array_temp,(3,natom,trotter))
2158      array_temp=zero
2159      do iimage=1,trotter
2160        iimagep=iimage+1;if(iimage==trotter) iimagep=1
2161        iimagem=iimage-1;if(iimage==1) iimagem=trotter
2162        do iatom=1,natom
2163          do ii=1,3
2164            array_temp(ii,iatom,iimage)=(dble(iimagem)*array(ii,iatom,iimagep)+array(ii,iatom,1))/dble(iimage)
2165          end do
2166        end do
2167      end do
2168      if (trotter>1) then
2169        do iimage=2,trotter
2170          do iatom=1,natom
2171            do ii=1,3
2172              array(ii,iatom,iimage)=array(ii,iatom,iimage)-array_temp(ii,iatom,iimage)
2173            end do
2174          end do
2175        end do
2176      end if
2177      ABI_FREE(array_temp)
2178 
2179 !  ---------------- From transformed to primitive coordinates ------------
2180    else if (ioption==-1) then
2181 
2182      ABI_MALLOC(array_temp,(3,natom,trotter))
2183      array_temp=zero
2184      do iimage=1,trotter
2185        do iatom=1,natom
2186          do ii=1,3
2187            array_temp(ii,iatom,iimage)=array(ii,iatom,1)
2188          end do
2189        end do
2190      end do
2191      if (trotter>1) then
2192        do iimage=2,trotter
2193          do iatom=1,natom
2194            do ii=1,3
2195              do iimage2=iimage,trotter
2196                array_temp(ii,iatom,iimage)=array_temp(ii,iatom,iimage) &
2197 &                     +array(ii,iatom,iimage2)*(dble(iimage-1))/(dble(iimage2-1))
2198              end do
2199            end do
2200          end do
2201        end do
2202      end if
2203      array(:,:,:)=array_temp(:,:,:)
2204      ABI_FREE(array_temp)
2205 
2206    end if ! ioption
2207 
2208  end if ! transform
2209 
2210 end subroutine pimd_coord_transform

m_pimd/pimd_destroy [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_destroy

FUNCTION

  Destroy the content of a datastructure of type pimd_type.
  Close open file(s) related to this datastructure.

INPUTS

OUTPUT

SIDE EFFECTS

  pimd_param=datastructure of type pimd_type.
            several parameters for PIMD.

SOURCE

275 subroutine pimd_destroy(pimd_param)
276 
277  implicit none
278 
279 !Arguments ------------------------------------
280 !scalars
281  type(pimd_type),intent(inout) :: pimd_param
282 !arrays
283 !Local variables-------------------------------
284 !scalars
285  integer :: ierr
286  character(len=100) :: msg
287 !arrays
288 
289 !************************************************************************
290 
291  ABI_SFREE(pimd_param%zeta_prev)
292  ABI_SFREE(pimd_param%zeta)
293  ABI_SFREE(pimd_param%zeta_next)
294  ABI_SFREE(pimd_param%dzeta)
295 
296  if (pimd_param%qtb_file_unit>0) then
297    if (is_open(pimd_param%qtb_file_unit)) then
298      ierr=close_unit(pimd_param%qtb_file_unit,msg)
299    end if
300  end if
301 
302  if (pimd_param%traj_unit>0) then
303    if (is_open(pimd_param%traj_unit)) then
304      ierr=close_unit(pimd_param%traj_unit,msg)
305    end if
306  end if
307 
308  call pimd_nullify(pimd_param)
309 
310 end subroutine pimd_destroy

m_pimd/pimd_diff_stress [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_diff_stress

FUNCTION

  Compute the difference between the stress tensor and the stress target

INPUTS

  stress_pimd(3,3,3)=stress tensor for PIMD
                     Last dimension (3) corresponds to 3 different pressure estimators
  stress_target(6)=stress target

OUTPUT

  pimd_diff_stress(3,3,3)=difference between stresses and stress target
                          First dimension (3) corresponds to 3 different pressure estimators

SIDE EFFECTS

SOURCE

1696 function pimd_diff_stress(stress_pimd,stress_target)
1697 
1698 !Arguments ------------------------------------
1699 !scalars
1700 !arrays
1701  real(dp),intent(in) :: stress_pimd(3,3,3),stress_target(6)
1702  real(dp) :: pimd_diff_stress(3,3)
1703 !Local variables-------------------------------
1704 !scalars
1705  integer :: ii,jj
1706 !arrays
1707  real(dp) :: stress_pimd2(3,3)
1708 
1709 !************************************************************************
1710 
1711 !Choice: the primitive estimator for pressure is chosen
1712  do ii=1,3
1713    do jj=1,3
1714      stress_pimd2(ii,jj)=stress_pimd(1,ii,jj)
1715    end do
1716  end do
1717 
1718 !+stress_target instead of - because it is translated from stress to pressure tensor
1719  pimd_diff_stress(1,1)=stress_pimd2(1,1)+stress_target(1)
1720  pimd_diff_stress(2,2)=stress_pimd2(2,2)+stress_target(2)
1721  pimd_diff_stress(3,3)=stress_pimd2(3,3)+stress_target(3)
1722  pimd_diff_stress(2,3)=stress_pimd2(2,3)+stress_target(4)
1723  pimd_diff_stress(3,2)=stress_pimd2(3,2)+stress_target(4)
1724  pimd_diff_stress(1,3)=stress_pimd2(1,3)+stress_target(5)
1725  pimd_diff_stress(3,1)=stress_pimd2(3,1)+stress_target(5)
1726  pimd_diff_stress(1,2)=stress_pimd2(1,2)+stress_target(6)
1727  pimd_diff_stress(2,1)=stress_pimd2(2,1)+stress_target(6)
1728 
1729 end function pimd_diff_stress

m_pimd/pimd_energies [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_energies

FUNCTION

  In the case od PIMD, compute the several contribution to total energy

INPUTS

  etotal_img(trotter)= energy (from DFT) for each cell
  forces(3,natom,trotter)=forces (from DFT) on atoms in each cell
  natom=number of atoms
  spring(natom)=spring constants in the primitive scheme
  trotter=Trotter number
  xcart(3,natom,trotter)=cartesian coordinates of atoms in each cell at t

OUTPUT

  eharm       =harmonic energy
  eharm_virial=harmonic energy from virial estimator
  epot        =potential energy

SIDE EFFECTS

SOURCE

1260 subroutine pimd_energies(eharm,eharm_virial,epot,etotal_img,forces,natom,spring,trotter,xcart)
1261 
1262 !Arguments ------------------------------------
1263 !scalars
1264  integer,intent(in) :: natom,trotter
1265  real(dp),intent(out) :: eharm,eharm_virial,epot
1266 !arrays
1267  real(dp),intent(in) :: etotal_img(trotter),forces(3,natom,trotter)
1268  real(dp),intent(in) :: xcart(3,natom,trotter)
1269  real(dp),intent(in) :: spring(natom)
1270 !Local variables-------------------------------
1271 !scalars
1272  integer :: iatom,ii,iimage,iimagep
1273 !arrays
1274  real(dp),allocatable :: centroid(:,:)
1275 
1276 !************************************************************************
1277 
1278 !Compute the centroid
1279  ABI_MALLOC(centroid,(3,natom))
1280  centroid=zero
1281  do iimage=1,trotter
1282    do iatom=1,natom
1283      do ii=1,3
1284        centroid(ii,iatom)=centroid(ii,iatom)+xcart(ii,iatom,iimage)
1285      end do
1286    end do
1287  end do
1288  centroid=centroid/dble(trotter)
1289 
1290 !Potential energy
1291  epot=sum(etotal_img(1:trotter))/dble(trotter)
1292 
1293 !Harmonic energy
1294  eharm=zero
1295  do iimage=1,trotter
1296    iimagep=iimage+1;if(iimage==trotter)iimagep=1
1297    do iatom=1,natom
1298      do ii=1,3
1299        eharm=eharm+half*spring(iatom)*((xcart(ii,iatom,iimagep)-xcart(ii,iatom,iimage))**2)
1300      end do
1301    end do
1302  end do
1303 
1304 !Harmonic energy from virial estimator
1305  eharm_virial=zero
1306  do iimage=1,trotter
1307    do iatom=1,natom
1308      do ii=1,3
1309        eharm_virial=eharm_virial-(xcart(ii,iatom,iimage)-centroid(ii,iatom)) &
1310 &              *forces(ii,iatom,iimage)
1311      end do
1312    end do
1313  end do
1314  eharm_virial=eharm_virial/dble(two*trotter)
1315 
1316  ABI_FREE(centroid)
1317 
1318 end subroutine pimd_energies

m_pimd/pimd_force_transform [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_force_transform

FUNCTION

  Apply a coordinate transformation on forces (defined for each atom in in cell) - for PIMD
  Possible choices for the transformation:
    0: no transformation
    1: normal mode tranformation
    2: staging transformation

INPUTS

  ioption=option given the direction of the transformation
     +1: from primitive coordinates to transformed coordinates
     -1: from transformed coordinates to primitive coordinates
  natom=number of atoms
  transform=coordinate transformation:
            0: no tranformation
            1: normal mode transformation
            2: staging transformation
  trotter=Trotter number

OUTPUT

SIDE EFFECTS

  forces(3,natom,trotter)=array containing forces

NOTES

  Back transformation (ioption=-1) not implemented !

SOURCE

2247 subroutine pimd_force_transform(forces,ioption,natom,transform,trotter)
2248 
2249 !Arguments ------------------------------------
2250 !scalars
2251  integer,intent(in) :: ioption,natom,transform,trotter
2252 !arrays
2253  real(dp),intent(inout) :: forces(3,natom,trotter)
2254 !Local variables-------------------------------
2255 !scalars
2256  integer :: iatom,ii,iimage,ll
2257  character(len=500) :: msg
2258 !arrays
2259  real(dp),allocatable :: forces_temp(:,:,:),nrm(:,:)
2260 
2261 !************************************************************************
2262 
2263  if (ioption==-1) then
2264    msg='Back transformation not implemented !'
2265    ABI_BUG(msg)
2266  end if
2267 
2268 !=== No transformation ===================================================
2269  if (transform==0) then
2270 
2271    return
2272 
2273 !=== Normal mode transformation ==========================================
2274  else if (transform==1) then
2275 
2276    !normal mode forces
2277    ABI_MALLOC(forces_temp,(3,natom,trotter))
2278    ABI_MALLOC(nrm,(trotter,trotter))
2279 
2280    do iimage=1,trotter
2281      nrm(iimage,1)=one
2282      nrm(iimage,trotter)=dble((-1)**(iimage-1))
2283    end do
2284 
2285    do ll=2,trotter/2
2286      do iimage=1,trotter
2287        nrm(iimage,2*ll-2)= two*cos(two_pi*dble((ll-1)* &
2288 &                         (iimage-1))/dble(trotter))
2289        nrm(iimage,2*ll-1)=-two*sin(two_pi*dble((ll-1)* &
2290 &                         (iimage-1))/dble(trotter))
2291      end do
2292    end do
2293 
2294    do ll=1,trotter
2295      forces_temp(:,:,ll)=zero
2296      do iimage=1,trotter
2297        forces_temp(:,:,ll)=forces_temp(:,:,ll)+nrm(iimage,ll)*forces(:,:,iimage)
2298      end do
2299    end do
2300 
2301    forces=forces_temp
2302 
2303    ABI_FREE(forces_temp)
2304    ABI_FREE(nrm)
2305 
2306 !=== Staging transformation ==============================================
2307  else if (transform==2) then
2308 
2309    !staging forces
2310    ABI_MALLOC(forces_temp,(3,natom,trotter))
2311    forces_temp=zero
2312    do iimage=1,trotter
2313      do iatom=1,natom
2314        do ii=1,3
2315          forces_temp(ii,iatom,1)=forces_temp(ii,iatom,1)+forces(ii,iatom,iimage)
2316        end do
2317      end do
2318    end do
2319    if (trotter>1) then
2320      do iimage=2,trotter
2321        do iatom=1,natom
2322          do ii=1,3
2323            forces_temp(ii,iatom,iimage)=forces(ii,iatom,iimage) &
2324 &               +forces_temp(ii,iatom,iimage-1)*(dble(iimage-2)/dble(iimage-1))
2325          end do
2326        end do
2327      end do
2328    end if
2329    forces=forces_temp
2330    ABI_FREE(forces_temp)
2331 
2332  end if ! transform
2333 
2334 end subroutine pimd_force_transform

m_pimd/pimd_forces [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_forces

FUNCTION

  Modify forces in order to take into account PIMD contribution

INPUTS

  natom=number of atoms
  spring(natom,spring_dim)=spring constants (spring_dim=1 or trotter)
  transform=coordinate transformation:
            0: no tranformation
            1: normal mode transformation
            2: staging transformation
  trotter=Trotter number
  xcart(3,natom,trotter)=cartesian coordinates of atoms in each cell

OUTPUT

SIDE EFFECTS

  forces(3,natom,trotter)=
    at input:  forces from electronic calculation
    at output: forces from electronic calculation + quantum spring contribution

SOURCE

1349 subroutine pimd_forces(forces,natom,spring,transform,trotter,xcart)
1350 
1351 !Arguments ------------------------------------
1352 !scalars
1353  integer,intent(in) :: natom,transform,trotter
1354 !arrays
1355  real(dp),intent(in) :: xcart(3,natom,trotter)
1356  real(dp),intent(in) :: spring(:,:)
1357  real(dp),intent(inout) :: forces(3,natom,trotter)
1358 !Local variables-------------------------------
1359 !scalars
1360  integer :: iatom,ii,iimage,iimagem,iimagep,ispring,natom_spring,nspring
1361  character(len=500) :: msg
1362 !arrays
1363 
1364 !************************************************************************
1365 
1366  natom_spring=size(spring,1);nspring=size(spring,2)
1367  if (natom/=natom_spring.or.(nspring/=1.and.nspring/=trotter)) then
1368    msg='Wrong dimensions for array spring !'
1369    ABI_BUG(msg)
1370  end if
1371 
1372  if (transform==0) then
1373    do iimage=1,trotter
1374      ispring=min(nspring,iimage)
1375      iimagep=iimage+1; iimagem=iimage-1
1376      if(iimage==trotter) iimagep=1
1377      if(iimage==1)       iimagem=trotter
1378      do iatom=1,natom
1379        do ii=1,3
1380          forces(ii,iatom,iimage)= &
1381 &               forces(ii,iatom,iimage)/dble(trotter) &
1382 &             - spring(iatom,ispring)*(two*xcart(ii,iatom,iimage)-xcart(ii,iatom,iimagem) &
1383 &                                                                -xcart(ii,iatom,iimagep))
1384        end do
1385      end do
1386    end do
1387 
1388  else
1389    do iimage=1,trotter
1390      ispring=min(nspring,iimage)
1391      do iatom=1,natom
1392        do ii=1,3
1393          forces(ii,iatom,iimage)= &
1394 &               forces(ii,iatom,iimage)/dble(trotter) &
1395 &             - spring(iatom,ispring)*xcart(ii,iatom,iimage)
1396        end do
1397      end do
1398    end do
1399 
1400  end if
1401 
1402 end subroutine pimd_forces

m_pimd/pimd_init [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_init

FUNCTION

  Initialize a datastructure of type pimd_type.
  Open file(s) related to this datastructure.

INPUTS

  dtset <type(dataset_type)>=all input variables in current dataset
  is_master= TRUE if I am the master process (proc 0)

OUTPUT

SIDE EFFECTS

  pimd_param=datastructure of type pimd_type.
             several parameters for Path-Integral MD.

SOURCE

135 subroutine pimd_init(dtset,pimd_param,is_master)
136 
137  implicit none
138 
139 !Arguments ------------------------------------
140 !scalars
141  logical,intent(in) :: is_master
142  type(dataset_type),target,intent(in) :: dtset
143  type(pimd_type),intent(inout) :: pimd_param
144 !Local variables-------------------------------
145 !scalars
146  integer :: ierr
147  character(len=200) :: msg
148 
149 !************************************************************************
150 
151  call pimd_nullify(pimd_param)
152 
153  if((dtset%imgmov==9).or.(dtset%imgmov==10).or.(dtset%imgmov==13))then
154    pimd_param%adpimd      = dtset%adpimd
155    pimd_param%constraint  = dtset%pimd_constraint
156    pimd_param%irandom     = dtset%irandom
157    pimd_param%nnos        = dtset%nnos
158    pimd_param%ntypat      = dtset%ntypat
159    pimd_param%optcell     = dtset%optcell
160    pimd_param%pitransform = dtset%pitransform
161    pimd_param%adpimd_gamma= dtset%adpimd_gamma
162    pimd_param%vis         = dtset%vis
163    pimd_param%bmass       = dtset%bmass
164    pimd_param%dtion       = dtset%dtion
165    pimd_param%friction    = dtset%friction
166    pimd_param%mdtemp      =>dtset%mdtemp
167    pimd_param%pimass      =>dtset%pimass
168    pimd_param%strtarget   =>dtset%strtarget
169    pimd_param%amu         =>dtset%amu_orig(:,1)
170    pimd_param%qmass       =>dtset%qmass
171    pimd_param%typat       =>dtset%typat
172    pimd_param%wtatcon     =>dtset%wtatcon
173    if(dtset%imgmov==10)then
174      pimd_param%use_qtb=1
175      if(is_master)then
176        call pimd_init_qtb(dtset,pimd_param%qtb_file_unit)
177      end if
178    end if
179    if(dtset%imgmov==13)then
180      ABI_MALLOC(pimd_param%zeta_prev,(3,dtset%natom,dtset%nimage,dtset%nnos))
181      ABI_MALLOC(pimd_param%zeta     ,(3,dtset%natom,dtset%nimage,dtset%nnos))
182      ABI_MALLOC(pimd_param%zeta_next,(3,dtset%natom,dtset%nimage,dtset%nnos))
183      ABI_MALLOC(pimd_param%dzeta    ,(3,dtset%natom,dtset%nimage,dtset%nnos))
184      pimd_param%zeta_prev=zero
185      pimd_param%zeta     =zero
186      pimd_param%zeta_next=zero
187      pimd_param%dzeta    =zero
188    end if
189    if(dtset%useria==37)then
190      ierr=open_file('pimd_traj.dat',msg,newunit=pimd_param%traj_unit,form='unformatted')
191      if (ierr/=0) then
192        ABI_ERROR(msg)
193      end if
194    end if
195  end if
196 
197 end subroutine pimd_init

m_pimd/pimd_init_qtb [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_init_qtb

FUNCTION

  Only relevant for PIMD + Quantum Thermal Bath (QTB);
  Initialize reading of PIQTB random force file.
  This routine should be called only by master proc.

INPUTS

  dtset <type(dataset_type)>=all input variables in current dataset

OUTPUT

  qtb_file_unit=if a PIQTB_force file exists, return its file unit.

SOURCE

332 subroutine pimd_init_qtb(dtset,qtb_file_unit)
333 
334  implicit none
335 
336 !Arguments ------------------------------------
337 !scalars
338  integer,intent(out) :: qtb_file_unit
339  type(dataset_type),target,intent(in) :: dtset
340 !Local variables-------------------------------
341 !scalars
342  integer :: ierr,ndof_qtb,ntimimage_qtb,nimage_qtb
343  real(dp) :: dtion_qtb,mdtemp_qtb
344  character(len=200) :: msg
345 
346 !************************************************************************
347 
348 !Try to open PIQTB random force file
349  ierr=open_file('piqtb_force',msg,newunit=qtb_file_unit,&
350 &               form='unformatted',status='old')
351 
352 !Read first line of the file
353  read(qtb_file_unit) dtion_qtb,ntimimage_qtb,mdtemp_qtb,nimage_qtb,ndof_qtb
354 
355 !Check consistency of the read parameters with ABINIT input file
356  if (abs(dtion_qtb-dtset%dtion)>tol6) then
357    msg='dtion read from piqtb_force file different from dtion in input file!'
358    ABI_ERROR(msg)
359  end if
360  if (abs(mdtemp_qtb-dtset%mdtemp(2))>tol6) then
361    msg='mdtemp read from piqtb_force file different from mdtemp(2) in input file!'
362    ABI_ERROR(msg)
363  end if
364  if (ntimimage_qtb<dtset%ntimimage) then
365    msg='ntimimage read from piqtb_force file smaller than ntimimage in input file!'
366    ABI_ERROR(msg)
367  end if
368  if (nimage_qtb/=dtset%nimage) then
369    msg='nimage read from piqtb_force file different from nimage in input file!'
370    ABI_ERROR(msg)
371  end if
372  if (ndof_qtb/=3*dtset%natom*dtset%nimage) then
373    msg='Nb of degrees of freedom read from piqtb_force not consistent with input file!'
374    ABI_ERROR(msg)
375  end if
376 
377 end subroutine pimd_init_qtb

m_pimd/pimd_initvel [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_initvel

FUNCTION

  Initialize velocities for PIMD with a gaussian distribution
  fixing the center of mass
  and eventually applying a constraint on atomic positions

INPUTS

  constraint=type of constraint to be applied
  mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
  natom=number of atoms
  temperature=temperature used to define velocities
  trotter=Trotter number
  wtatcon(3,natom)=weights for atomic constraints

OUTPUT

  vel(3,natom,trotter)=velocities of atoms in in each cell

SIDE EFFECTS

  iseed=seed for random number generator

SOURCE

838 subroutine pimd_initvel(iseed,mass,natom,temperature,trotter,vel,constraint,wtatcon)
839 
840 !Arguments ------------------------------------
841 !scalars
842  integer,intent(in) :: constraint,natom,trotter
843  integer,intent(inout) :: iseed
844  real(dp),intent(in) :: temperature
845 !arrays
846  real(dp),intent(in) :: mass(:,:),wtatcon(3,natom)
847  real(dp),intent(out) :: vel(3,natom,trotter)
848 !Local variables-------------------------------
849 !scalars
850  integer :: iatom,ii,iimage,imass,natom_mass,nmass
851  real(dp) :: mtot,rescale_vel
852  character(len=500) :: msg
853 !arrays
854  real(dp) :: mvini(3)
855 
856 !************************************************************************
857 
858  natom_mass=size(mass,1);nmass=size(mass,2)
859  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
860    msg='Wrong dimensions for array mass !'
861    ABI_BUG(msg)
862  end if
863 
864 !Compute total mass (of non-constrained atoms)
865  if (constraint==0) then
866    mtot=sum(mass(1:natom,1:nmass))
867  else
868    mtot=zero
869    do iatom=1,natom
870      if (all(abs(wtatcon(:,iatom))<tol8)) mtot=mtot+sum(mass(iatom,1:nmass))
871    end do
872  end if
873  if (nmass==1) mtot=mtot*dble(trotter)
874 
875 !Initialize randomly the velocities
876  do iimage=1,trotter
877    imass=min(nmass,iimage)
878    do iatom=1,natom
879      do ii=1,3
880        vel(ii,iatom,iimage)=sqrt(kb_HaK*temperature/mass(iatom,imass))*cos(two_pi*uniformrandom(iseed))
881        vel(ii,iatom,iimage)=vel(ii,iatom,iimage)*sqrt(-two*log(uniformrandom(iseed)))
882      end do
883    end do
884  end do
885 
886 !Cancel velocities of constrained atoms
887  if (constraint/=0) then
888    do iimage=1,trotter
889      do iatom=1,natom
890        if (any(abs(wtatcon(:,iatom))>=tol8)) vel(:,iatom,iimage)=zero
891      end do
892    end do
893  end if
894 
895 !Make sure that the (sum of m_i v_i) at step zero is zero
896  mvini=zero
897  do iimage=1,trotter
898    imass=min(nmass,iimage)
899    do iatom=1,natom
900      do ii=1,3
901       mvini(ii)=mvini(ii)+mass(iatom,imass)*vel(ii,iatom,iimage)
902      end do
903    end do
904  end do
905  if (constraint==0) then
906    do iimage=1,trotter
907      do iatom=1,natom
908        do ii=1,3
909          vel(ii,iatom,iimage)=vel(ii,iatom,iimage)-(mvini(ii)/mtot)
910        end do
911      end do
912    end do
913  else
914    do iimage=1,trotter
915      do iatom=1,natom
916        if (all(abs(wtatcon(:,iatom))<tol8)) vel(:,iatom,iimage)=vel(:,iatom,iimage)-(mvini(:)/mtot)
917      end do
918    end do
919  end if
920 
921 !Now rescale the velocities to give the exact temperature
922  rescale_vel=sqrt(temperature/pimd_temperature(mass,vel))
923  vel(:,:,:)=vel(:,:,:)*rescale_vel
924 
925 end subroutine pimd_initvel

m_pimd/pimd_is_restart [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_is_restart

FUNCTION

  Determine whether this is a PIMD restart or not:
  test on value of velocities and corresponding temperature

INPUTS

  mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
  vel(3,natom,nimage)=velocities for each image of the cell

OUTPUT

  pimd_is_restart=1 if temperature is not zero

SIDE EFFECTS

SOURCE

445 function pimd_is_restart(mass,vel,vel_cell)
446 
447 !Arguments ------------------------------------
448 !scalars
449  integer :: pimd_is_restart
450 !arrays
451  real(dp),intent(in) :: mass(:,:),vel(:,:,:)
452  real(dp),intent(in),optional :: vel_cell(:,:)
453 !Local variables-------------------------------
454 !scalars
455  real(dp),parameter :: zero_temp=tol7
456 !arrays
457 
458 !************************************************************************
459 
460  pimd_is_restart=0
461  if (pimd_temperature(mass,vel)>zero_temp) pimd_is_restart=1
462  if (present(vel_cell)) then
463    if (maxval(vel_cell)>zero_temp) pimd_is_restart=pimd_is_restart+10
464  end if
465 
466 end function pimd_is_restart

m_pimd/pimd_langevin_forces [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_langevin_forces

FUNCTION

  Compute Langevin contribution to PIMD forces

INPUTS

  alea(3,natom,trotter)=set of random numbers
  forces(3,natom,trotter)=forces without Langevin contribution
  friction=friction factor
  langev(natom,mass_dim)=Langevin factors (mass_dim=1 or trotter)
  mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
  natom=number of atoms
  trotter=Trotter number
  vel(3,natom,trotter)=velocities of atoms in each cell

OUTPUT

  forces_langevin(3,natom,trotter)=forces including Langevin contribution

SIDE EFFECTS

SOURCE

1431 subroutine pimd_langevin_forces(alea,forces,forces_langevin,friction,&
1432 &                               langev,mass,natom,trotter,vel)
1433 
1434 !Arguments ------------------------------------
1435 !scalars
1436  integer,intent(in) :: natom,trotter
1437  real(dp),intent(in) :: friction
1438 !arrays
1439  real(dp),intent(in) :: alea(3,natom,trotter),forces(3,natom,trotter)
1440  real(dp),intent(in) :: langev(:,:),mass(:,:),vel(3,natom,trotter)
1441  real(dp),intent(out) :: forces_langevin(3,natom,trotter)
1442 !Local variables-------------------------------
1443 !scalars
1444  integer :: iatom,ii,iimage,imass,natom_mass,nmass
1445  character(len=500) :: msg
1446 !arrays
1447 
1448 !************************************************************************
1449 
1450  natom_mass=size(mass,1);nmass=size(mass,2)
1451  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
1452    msg='Wrong dimensions for array mass !'
1453    ABI_BUG(msg)
1454  end if
1455 
1456  do iimage=1,trotter
1457    imass=min(nmass,iimage)
1458    do iatom=1,natom
1459      do ii=1,3
1460        forces_langevin(ii,iatom,iimage)=forces(ii,iatom,iimage) &
1461 &                    + langev(iatom,imass)*alea(ii,iatom,iimage) &
1462 &                    - friction*mass(iatom,imass)*vel(ii,iatom,iimage)
1463      end do
1464    end do
1465  end do
1466 
1467 end subroutine pimd_langevin_forces

m_pimd/pimd_langevin_random [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_langevin_random

FUNCTION

  Generate a set of random numbers to be used for PIMD Langevin algorithm

INPUTS

  irandom=option for random number generator:
          1:uniform random routine provided within Abinit package
          2:Fortran 90 random number generator
          3:ZBQLU01 non deterministic random number generator
  langev(natom,mass_dim)=Langevin factors (mass_dim=1 or trotter)
  mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
  natom=number of atoms
  trotter=Trotter number
  zeroforce=flag; if 1 keep sum of forces equal to zero

OUTPUT

  alea(3,natom,trotter)=set of random numbers

SIDE EFFECTS

  iseed=seed for random number generator (used only if irandom=1)

SOURCE

 956 subroutine pimd_langevin_random(alea,irandom,iseed,langev,mass,natom,trotter,zeroforce)
 957 
 958 !Arguments ------------------------------------
 959 !scalars
 960  integer,intent(in) :: irandom,natom,trotter,zeroforce
 961  integer,intent(inout) :: iseed
 962 !arrays
 963  real(dp),intent(in) :: langev(:,:),mass(:,:)
 964  real(dp),intent(out) :: alea(3,natom,trotter)
 965 !Local variables-------------------------------
 966 !scalars
 967  integer :: iatom,ii,iimage,imass,nmass,natom_mass
 968  real(dp) :: mtot,r1,r2
 969  character(len=500) :: msg
 970 !arrays
 971  real(dp) :: total(3)
 972 
 973 !************************************************************************
 974 
 975  natom_mass=size(mass,1);nmass=size(mass,2)
 976  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
 977    msg='Wrong dimensions for array mass !'
 978    ABI_BUG(msg)
 979  end if
 980 
 981  mtot=sum(mass(1:natom,1:nmass))
 982  if (nmass==1) mtot=mtot*trotter
 983 
 984 !Draw random numbers
 985  do iimage=1,trotter
 986    do iatom=1,natom
 987      do ii=1,3
 988        select case(irandom)
 989        case(1)
 990          r1=uniformrandom(iseed)
 991          r2=uniformrandom(iseed)
 992        case(2)
 993          call random_number(r1)
 994          call random_number(r2)
 995        case(3)
 996          r1=ZBQLU01(zero)
 997          r2=ZBQLU01(zero)
 998        end select
 999        alea(ii,iatom,iimage)= cos(two_pi*r1)*sqrt(-log(r2)*two)
1000      end do
1001    end do
1002  end do
1003 
1004 !Make sure that the sum of random forces is zero
1005  if(zeroforce==1)then
1006    total=zero
1007    mtot=zero
1008    do iimage=1,trotter
1009      imass=min(nmass,iimage)
1010      do iatom=1,natom
1011        mtot=mtot+mass(iatom,imass)
1012      end do
1013    end do
1014    do iimage=1,trotter
1015      imass=min(nmass,iimage)
1016      do iatom=1,natom
1017        do ii=1,3
1018          total(ii)=total(ii)+langev(iatom,imass)*alea(ii,iatom,iimage)
1019        end do
1020      end do
1021    end do
1022    do iimage=1,trotter
1023      imass=min(nmass,iimage)
1024      do iatom=1,natom
1025        do ii=1,3
1026          alea(ii,iatom,iimage)= alea(ii,iatom,iimage)- &
1027 &        (total(ii)*mass(iatom,imass))/(langev(iatom,imass)*mtot)
1028        end do
1029      end do
1030    end do
1031 !  now random forces have been rescaled so that their sum is zero
1032  end if
1033 
1034 end subroutine pimd_langevin_random

m_pimd/pimd_langevin_random_bar [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_langevin_random_bar

FUNCTION

  Generate a set of random numbers to be used for the barostat of PIMD Langevin algorithm

INPUTS

  irandom=option for random number generator:
          1:uniform random routine provided within Abinit package
          2:Fortran 90 random number generator
          3:ZBQLU01 non deterministic random number generator

OUTPUT

  alea_bar(3,3)=set of random numbers

SIDE EFFECTS

  iseed=seed for random number generator (used only if irandom=1)

SOURCE

1148 subroutine pimd_langevin_random_bar(alea_bar,irandom,iseed)
1149 
1150 !Arguments ------------------------------------
1151 !scalars
1152  integer,intent(in) :: irandom
1153  integer,intent(inout) :: iseed
1154 !arrays
1155  real(dp),intent(out) :: alea_bar(3,3)
1156 !Local variables-------------------------------
1157 !scalars
1158  integer :: ii,jj
1159  real(dp) :: r1,r2
1160 !arrays
1161 
1162 !************************************************************************
1163 
1164 !Draw random numbers
1165  do ii=1,3
1166    do jj=1,3
1167      select case(irandom)
1168      case(1)
1169        r1=uniformrandom(iseed)
1170        r2=uniformrandom(iseed)
1171      case(2)
1172        call random_number(r1)
1173        call random_number(r2)
1174      case(3)
1175        r1=ZBQLU01(zero)
1176        r2=ZBQLU01(zero)
1177      end select
1178      alea_bar(ii,jj)= cos(two_pi*r1)*sqrt(-log(r2)*two)
1179    end do
1180  end do
1181 
1182 !Symmetrize
1183  alea_bar(1,2)=half*(alea_bar(1,2)+alea_bar(2,1))
1184  alea_bar(1,3)=half*(alea_bar(1,3)+alea_bar(3,1))
1185  alea_bar(2,3)=half*(alea_bar(2,3)+alea_bar(3,2))
1186  alea_bar(2,1)=alea_bar(1,2)
1187  alea_bar(3,1)=alea_bar(1,3)
1188  alea_bar(3,2)=alea_bar(2,3)
1189 
1190 end subroutine pimd_langevin_random_bar

m_pimd/pimd_langevin_random_init [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_langevin_random_init

FUNCTION

  Initialize random number generator to be used for PIMD Langevin algorithm

INPUTS

  irandom=option for random number generator:
          1:uniform random routine provided within Abinit package
          2:Fortran 90 random number generator
          3:ZBQLU01 non deterministic random number generator

OUTPUT

SIDE EFFECTS

  iseed=seed for random number generator (used only if irandom=1)

SOURCE

1215 subroutine pimd_langevin_random_init(irandom,iseed)
1216 
1217 !Arguments ------------------------------------
1218 !scalars
1219  integer,intent(in) :: irandom
1220  integer,intent(inout) :: iseed
1221 
1222 !************************************************************************
1223 
1224  if (irandom==3) then
1225    call ZBQLINI(0)
1226  end if
1227 
1228 !Fake statement
1229  return;if (.false.) iseed=zero
1230 
1231 end subroutine pimd_langevin_random_init

m_pimd/pimd_langevin_random_qtb [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_langevin_random_qtb

FUNCTION

  Read a set of random forces (atm units) to be used for PIMD QTB algorithm

INPUTS

  langev(natom,mass_dim)=Langevin factors (mass_dim=1 or trotter)
  mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
  natom=number of atoms
  qtb_file_unit= ramdom forces file unit
  trotter=Trotter number
  zeroforce=flag; if 1 keep sum of forces equal to zero

OUTPUT

  alea(3,natom,trotter)=set of random forces

SOURCE

1059 subroutine pimd_langevin_random_qtb(alea,langev,mass,natom,qtb_file_unit,trotter,zeroforce)
1060 
1061 !Arguments ------------------------------------
1062 !scalars
1063  integer,intent(in) :: natom,qtb_file_unit,trotter,zeroforce
1064 !arrays
1065  real(dp),intent(in) :: langev(:,:),mass(:,:)
1066  real(dp),intent(out) :: alea(3,natom,trotter)
1067 !Local variables-------------------------------
1068 !scalars
1069  integer :: iatom,ii,iimage,imass,nmass,natom_mass
1070  real(dp) :: mtot
1071  character(len=500) :: msg
1072 !arrays
1073  real(sp) :: alea_sp(3,natom,trotter)
1074  real(dp) :: total(3)
1075 
1076 !************************************************************************
1077 
1078  if (qtb_file_unit<0) then
1079    msg='QTB forces file unit should be positive!'
1080    ABI_BUG(msg)
1081  end if
1082  natom_mass=size(mass,1);nmass=size(mass,2)
1083  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
1084    msg='Wrong dimensions for array mass !'
1085    ABI_BUG(msg)
1086  end if
1087 
1088 !Read QTB random forces
1089  read(qtb_file_unit) alea_sp(1:3,1:natom,1:trotter)
1090  alea(:,:,:)=dble(alea_sp(:,:,:))
1091 
1092 !Make sure that the sum of random forces is zero
1093  if(zeroforce==1)then
1094    total=zero
1095    mtot=zero
1096    do iimage=1,trotter
1097      imass=min(nmass,iimage)
1098      do iatom=1,natom
1099        mtot=mtot+mass(iatom,imass)
1100      end do
1101    end do
1102    do iimage=1,trotter
1103      imass=min(nmass,iimage)
1104      do iatom=1,natom
1105        do ii=1,3
1106          total(ii)=total(ii)+langev(iatom,imass)*alea(ii,iatom,iimage)
1107        end do
1108      end do
1109    end do
1110    do iimage=1,trotter
1111      imass=min(nmass,iimage)
1112      do iatom=1,natom
1113        do ii=1,3
1114          alea(ii,iatom,iimage)= alea(ii,iatom,iimage)- &
1115 &        (total(ii)*mass(iatom,imass))/(langev(iatom,imass)*mtot)
1116        end do
1117      end do
1118    end do
1119 !  now random forces have been rescaled so that their sum is zero
1120  end if
1121 
1122 end subroutine pimd_langevin_random_qtb

m_pimd/pimd_mass_spring [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_mass_spring

FUNCTION

  Compute masses and spring constants for PIMD. Eventually apply a coordinate transformation.
  Possible choices for the transformation:
    0: no transformation
    1: normal mode tranformation
    2: staging transformation

INPUTS

  inertmass(natom)=fictitious masses of atoms
  kt=kT constant
  natom=number of atoms
  quantummass(natom)=true masses of atoms
  transform=coordinate transformation:
            0: no tranformation
            1: normal mode transformation
            2: staging transformation
  trotter=Trotter number

OUTPUT

SIDE EFFECTS

  mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
  spring(natom,mass_dim)=spring constants of atoms (mass_dim=1 or trotter)

NOTES

  Back transformation (ioption=-1) not implemented !

SOURCE

2493 subroutine pimd_mass_spring(inertmass,kt,mass,natom,quantummass,spring,transform,trotter)
2494 
2495 !Arguments ------------------------------------
2496 !scalars
2497  integer,intent(in) :: natom,transform,trotter
2498  real(dp),intent(in) :: kt
2499 !arrays
2500  real(dp),intent(in) :: inertmass(natom),quantummass(natom)
2501  real(dp),intent(out) :: mass(:,:),spring(:,:)
2502 !Local variables-------------------------------
2503 !scalars
2504  integer :: iimage,kk,natom_mass,natom_spring,nmass,nspring
2505  real(dp) :: gammasquare
2506  character(len=500) :: msg
2507 !arrays
2508  real(dp),allocatable :: mass_temp(:,:),lambda(:)
2509 
2510 !************************************************************************
2511 
2512  natom_mass  =size(mass  ,1);nmass  =size(mass  ,2)
2513  natom_spring=size(spring,1);nspring=size(spring,2)
2514  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
2515    msg='Wrong dimensions for array mass !'
2516    ABI_BUG(msg)
2517  end if
2518  if (natom/=natom_spring.or.(nspring/=1.and.nspring/=trotter)) then
2519    msg='Wrong dimensions for array spring !'
2520    ABI_BUG(msg)
2521  end if
2522 
2523 !=== No transformation ===================================================
2524  if (transform==0) then
2525    !2nd dim of mass and spring = 1
2526    mass(1:natom,1)=inertmass(1:natom)
2527    spring(1:natom,1)=quantummass(1:natom)*dble(trotter)*kt*kt
2528 
2529 !=== Normal mode transformation ==========================================
2530  else if (transform==1) then
2531 
2532    ABI_MALLOC(lambda,(trotter))
2533    lambda(1)=zero; lambda(trotter)=four*dble(trotter)
2534    do kk=2,trotter/2
2535      lambda(2*kk-2)=four*dble(trotter)* &
2536 &                  (one-cos(two_pi*dble(kk-1)/dble(trotter)))
2537      lambda(2*kk-1)=four*dble(trotter)* &
2538 &                  (one-cos(two_pi*dble(kk-1)/dble(trotter)))
2539    end do
2540 
2541    !normal mode masses
2542    do iimage=1,trotter
2543      mass(:,iimage)=quantummass(:)*lambda(iimage)
2544    end do
2545 
2546    do iimage=1,trotter
2547      spring(:,iimage)=mass(:,iimage)*dble(trotter)*kt*kt
2548    end do
2549 
2550    !fictitious masses
2551    mass(:,1)=inertmass(:)
2552 
2553    !from 2 to P not changed except if adiabatic PIMD
2554    !see Hone et al, JCP 124, 154103 (2006)
2555    gammasquare=one  !adiabaticity parameter
2556    do iimage=2,trotter
2557      mass(:,iimage)=mass(:,iimage)/gammasquare
2558    end do
2559 
2560    ABI_FREE(lambda)
2561 
2562 !=== Staging transformation ==============================================
2563  else if (transform==2) then
2564 
2565    !Fictitious masses
2566    mass(1:natom,1)=inertmass(1:natom)
2567    if (nmass>1) then
2568      do iimage=2,trotter
2569        mass(1:natom,iimage)=inertmass(1:natom)*dble(iimage)/dble(iimage-1)
2570      end do
2571    end if
2572 
2573    !Staging masses (mass_temp)
2574    ABI_MALLOC(mass_temp,(natom,trotter))
2575    mass_temp(1:natom,1)=zero
2576    if (nmass>1) then
2577      do iimage=2,trotter
2578        mass_temp(1:natom,iimage)=quantummass(1:natom)*dble(iimage)/dble(iimage-1)
2579      end do
2580    end if
2581 
2582    spring(1:natom,1)=mass_temp(1:natom,1)*dble(trotter)*kt*kt
2583    if (nspring>1) then
2584      do iimage=2,trotter
2585        spring(1:natom,iimage)=mass_temp(1:natom,iimage)*dble(trotter)*kt*kt
2586      end do
2587    end if
2588    ABI_FREE(mass_temp)
2589 
2590  end if
2591 
2592 end subroutine pimd_mass_spring

m_pimd/pimd_noseehoover_forces [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_nosehoover_forces

FUNCTION

  Compute Nose-Hoover contribution to PIMD forces
  by adding friction force of thermostat number one

INPUTS

  dzeta(3,natom,trotter,nnos)=variables of thermostats, in (atomic time unit)^(-1)
                              used only when a coordinate transformation is applied (transfom/=0)
  forces(3,natom,trotter)=forces without Nose-Hoover contribution
  mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
  natom=number of atoms
  nnos=number of thermostats
  trotter=Trotter number
  vel(3,natom,trotter)=velocities of atoms in each cell

OUTPUT

  forces_nosehoover(3,natom,trotter)=forces including thermostat contribution

SIDE EFFECTS

SOURCE

1497 subroutine pimd_nosehoover_forces(dzeta,forces,forces_nosehoover,mass,natom,&
1498 &                                 nnos,trotter,vel)
1499 
1500 !Arguments ------------------------------------
1501 !scalars
1502  integer,intent(in) :: natom,nnos,trotter
1503 !arrays
1504  real(dp),intent(in) :: dzeta(3,natom,trotter,nnos),forces(3,natom,trotter)
1505  real(dp),intent(in) :: vel(3,natom,trotter)
1506  real(dp),intent(in) :: mass(:,:)
1507  real(dp),intent(out) :: forces_nosehoover(3,natom,trotter)
1508 !Local variables-------------------------------
1509 !scalars
1510  integer :: iatom,ii,iimage,imass,natom_mass,nmass
1511  character(len=500) :: msg
1512 !arrays
1513 
1514 !************************************************************************
1515 
1516  natom_mass=size(mass,1);nmass=size(mass,2)
1517  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
1518    msg='Wrong dimensions for array mass !'
1519    ABI_BUG(msg)
1520  end if
1521 
1522   do iimage=1,trotter
1523     imass=min(nmass,iimage)
1524     do iatom=1,natom
1525       do ii=1,3
1526         forces_nosehoover(ii,iatom,iimage)=forces(ii,iatom,iimage) &
1527 &           - mass(iatom,imass)*dzeta(ii,iatom,iimage,1)*vel(ii,iatom,iimage)
1528       end do
1529     end do
1530   end do
1531 
1532 end subroutine pimd_nosehoover_forces

m_pimd/pimd_nosehoover_propagate [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_nosehoover_propagate

FUNCTION

  Propagate thermostat variables (Nose-Hoover algorithm) - for PIMD

INPUTS

  dtion=time step
  dzeta(3,natom,trotter,nnos)=time derivative of zeta at time t
  itimimage=index of time step
  mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
  natom=number of atoms
  nnos=number of thermostats
  qmass(nnos)=masses of thermostats
  temperature=temperature
  transform=coordinate transformation:
            0: no tranformation
            1: normal mode transformation
            2: staging transformation
  trotter=Trotter number
  vel(3,natom,trotter)=velocities of atoms in each cell
  zeta(3,natom,trotter,nnos)=variables of thermostats, in (atomic time unit)^(-1)
             used only when no coordinate transformation is applied (transfom==0)
             at time t
  zeta_prev(3,natom,trotter,nnos)=previous value of zeta (t-dt)

OUTPUT

  zeta_next(3,natom,trotter,nnos)=next value of zeta (t+dt)

SIDE EFFECTS

SOURCE

1896 subroutine pimd_nosehoover_propagate(dtion,dzeta,mass,natom,nnos,qmass,temperature,&
1897 &                                    trotter,vel,zeta,zeta_next,zeta_prev,itimimage,transform)
1898 
1899 !Arguments ------------------------------------
1900 !scalars
1901  integer,intent(in) :: itimimage,natom,nnos,transform,trotter
1902  real(dp),intent(in) :: dtion,temperature
1903 !arrays
1904  real(dp),intent(in) :: qmass(nnos),vel(3,natom,trotter)
1905  real(dp),intent(in) :: mass(:,:)
1906  real(dp),intent(in) :: dzeta(3,natom,trotter,nnos),zeta_prev(3,natom,trotter,nnos)
1907  real(dp),intent(in) :: zeta(3,natom,trotter,nnos)
1908  real(dp),intent(out) :: zeta_next(3,natom,trotter,nnos)
1909 !Local variables-------------------------------
1910 !scalars
1911  integer :: iatom,ii,iimage,inos,natom_mass,nmass
1912  real(dp) :: kt
1913  character(len=500) :: msg
1914 !arrays
1915  real(dp),allocatable :: thermforces(:,:,:,:)
1916 
1917 !************************************************************************
1918 
1919  natom_mass=size(mass,1);nmass=size(mass,2)
1920  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
1921    msg='Wrong dimensions for array mass !'
1922    ABI_BUG(msg)
1923  end if
1924  if (nnos<3) then
1925    msg='Not available for nnos<3 !'
1926    ABI_BUG(msg)
1927  end if
1928 
1929  kt=temperature*kb_HaK
1930  ABI_MALLOC(thermforces,(3,natom,trotter,nnos))
1931 
1932 !forces on the thermostats
1933  do inos=1,nnos
1934    if(inos==1)then
1935      do iimage=1,trotter
1936        do iatom=1,natom
1937          do ii=1,3
1938            select case(transform)
1939            case(0)
1940              thermforces(ii,iatom,iimage,1)=mass(iatom,1)*(vel(ii,iatom,iimage)**2)-kt
1941            case(1)
1942              thermforces(ii,iatom,iimage,1)=mass(iatom,iimage)*(vel(ii,iatom,iimage)**2)-kt
1943            case(2)
1944              thermforces(ii,iatom,iimage,1)=mass(iatom,iimage)*(vel(ii,iatom,iimage)**2)-kt
1945            end select
1946          end do
1947        end do
1948      end do
1949    elseif(inos==(nnos-1))then
1950      do iimage=1,trotter
1951        do iatom=1,natom
1952          do ii=1,3
1953            thermforces(ii,iatom,iimage,inos)=&
1954 &          (qmass(nnos-2)*(dzeta(ii,iatom,iimage,nnos-2)**2)-kt) + &
1955 &          (qmass(nnos  )*(dzeta(ii,iatom,iimage,nnos)  **2)-kt)
1956          end do
1957        end do
1958      end do
1959    elseif(inos==nnos)then
1960      do iimage=1,trotter
1961        do iatom=1,natom
1962          do ii=1,3
1963            thermforces(ii,iatom,iimage,inos)=&
1964 &          qmass(nnos-1)*(dzeta(ii,iatom,iimage,nnos-1)**2)-kt
1965          end do
1966        end do
1967      end do
1968    else
1969      do iimage=1,trotter
1970        do iatom=1,natom
1971          do ii=1,3
1972            thermforces(ii,iatom,iimage,inos)=&
1973 &          qmass(inos-1)*(dzeta(ii,iatom,iimage,inos-1)**2)-kt
1974          end do
1975        end do
1976      end do
1977    end if
1978  end do
1979 
1980  select case(itimimage)
1981  case(1) !taylor
1982 
1983  do inos=1,nnos
1984    if(inos==1)then
1985      zeta_next(:,:,:,1)=zeta(:,:,:,1)+ dzeta(:,:,:,1)*dtion + &
1986 &       (thermforces(:,:,:,1)-qmass(1)*dzeta(:,:,:,1)*dzeta(:,:,:,2))* &
1987 &       dtion*dtion/(two*qmass(1))
1988    elseif(inos==(nnos-1))then
1989      zeta_next(:,:,:,inos)=zeta(:,:,:,inos)+ dzeta(:,:,:,inos)*dtion + &
1990 &       (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,inos)*dzeta(:,:,:,nnos))* &
1991 &       dtion*dtion/(two*qmass(inos))
1992    elseif(inos==nnos)then
1993      zeta_next(:,:,:,inos)=zeta(:,:,:,inos)+ dzeta(:,:,:,inos)*dtion + &
1994 &       (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,nnos-1)*dzeta(:,:,:,nnos))* &
1995 &       dtion*dtion/(two*qmass(inos))
1996    else
1997      zeta_next(:,:,:,inos)=zeta(:,:,:,inos)+ dzeta(:,:,:,inos)*dtion + &
1998 &       (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,inos)*dzeta(:,:,:,inos+1))* &
1999 &       dtion*dtion/(two*qmass(inos))
2000    end if
2001  end do
2002 
2003  case default !verlet
2004 
2005  do inos=1,nnos
2006    if(inos==1)then
2007      zeta_next(:,:,:,1)=two*zeta(:,:,:,1) - zeta_prev(:,:,:,1) + &
2008 &       (thermforces(:,:,:,1)-qmass(1)*dzeta(:,:,:,1)*dzeta(:,:,:,2))* &
2009 &       dtion*dtion/qmass(1)
2010    elseif(inos==(nnos-1))then
2011      zeta_next(:,:,:,inos)=two*zeta(:,:,:,inos) - zeta_prev(:,:,:,inos) + &
2012 &       (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,inos)*dzeta(:,:,:,nnos))* &
2013 &       dtion*dtion/qmass(inos)
2014    elseif(inos==nnos)then
2015      zeta_next(:,:,:,inos)=two*zeta(:,:,:,inos) - zeta_prev(:,:,:,inos) + &
2016 &       (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,nnos-1)*dzeta(:,:,:,nnos))* &
2017 &       dtion*dtion/qmass(inos)
2018    else
2019      zeta_next(:,:,:,inos)=two*zeta(:,:,:,inos) - zeta_prev(:,:,:,inos) + &
2020 &       (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,inos)*dzeta(:,:,:,inos+1))* &
2021 &       dtion*dtion/qmass(inos)
2022    end if
2023  end do
2024 
2025  end select
2026 
2027  ABI_FREE(thermforces)
2028 
2029 end subroutine pimd_nosehoover_propagate

m_pimd/pimd_nullify [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_nullify

FUNCTION

  Nullify the content of a datastructure of type pimd_type.

INPUTS

OUTPUT

SIDE EFFECTS

  pimd_param=datastructure of type pimd_type.
             several parameters for Path-Integral MD.

SOURCE

219 subroutine pimd_nullify(pimd_param)
220 
221  implicit none
222 
223 !Arguments ------------------------------------
224 !scalars
225  type(pimd_type),intent(inout) :: pimd_param
226 
227 !************************************************************************
228 
229  pimd_param%adpimd       =  0
230  pimd_param%constraint   =  0
231  pimd_param%irandom      = -1
232  pimd_param%nnos         = -1
233  pimd_param%ntypat       = -1
234  pimd_param%optcell      = -1
235  pimd_param%pitransform  = -1
236  pimd_param%qtb_file_unit= -1
237  pimd_param%traj_unit    = -1
238  pimd_param%use_qtb      =  0
239  pimd_param%adpimd_gamma = one
240  pimd_param%vis          = zero
241  pimd_param%bmass        = zero
242  pimd_param%dtion        = zero
243  pimd_param%friction     = zero
244  nullify(pimd_param%mdtemp)
245  nullify(pimd_param%pimass)
246  nullify(pimd_param%strtarget)
247  nullify(pimd_param%amu)
248  nullify(pimd_param%qmass)
249  nullify(pimd_param%typat)
250  nullify(pimd_param%wtatcon)
251 
252 end subroutine pimd_nullify

m_pimd/pimd_predict_taylor [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_predict_taylor

FUNCTION

  Predict new atomic positions using a Taylor algorithm (first time step) - for PIMD

INPUTS

  dtion=time step
  forces(3,natom,trotter)=PIMD forces on atoms in each cell
  mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
  natom=number of atoms
  trotter=Trotter number
  vel(3,natom,trotter)=velocities of atoms in each cell
  xcart(3,natom,trotter)=cartesian coordinates of atoms in each cell at t

OUTPUT

  xcart_next(3,natom,trotter)=cartesian coordinates of atoms in each cell at t+dt

SIDE EFFECTS

SOURCE

1757 subroutine pimd_predict_taylor(dtion,forces,mass,natom,trotter,vel,xcart,xcart_next)
1758 
1759 !Arguments ------------------------------------
1760 !scalars
1761  integer,intent(in) :: natom,trotter
1762  real(dp),intent(in) :: dtion
1763 !arrays
1764  real(dp),intent(in) :: forces(3,natom,trotter),vel(3,natom,trotter),xcart(3,natom,trotter)
1765  real(dp),intent(in) :: mass(:,:)
1766  real(dp),intent(out) :: xcart_next(3,natom,trotter)
1767 !Local variables-------------------------------
1768 !scalars
1769  integer :: iatom,ii,iimage,imass,natom_mass,nmass
1770  character(len=500) :: msg
1771 !arrays
1772 
1773 !************************************************************************
1774 
1775  natom_mass=size(mass,1);nmass=size(mass,2)
1776  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
1777    msg='Wrong dimensions for array mass !'
1778    ABI_BUG(msg)
1779  end if
1780 
1781  do iimage=1,trotter
1782    imass=min(nmass,iimage)
1783    do iatom=1,natom
1784      do ii=1,3
1785        xcart_next(ii,iatom,iimage)=xcart(ii,iatom,iimage) &
1786 &             + half*dtion*dtion*forces(ii,iatom,iimage)/mass(iatom,imass) &
1787 &             + dtion*vel(ii,iatom,iimage)
1788      end do
1789    end do
1790  end do
1791 
1792 end subroutine pimd_predict_taylor

m_pimd/pimd_predict_verlet [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_predict_verlet

FUNCTION

  Predict new atomic positions using a Verlet algorithm - for PIMD

INPUTS

  dtion=time step
  forces(3,natom,trotter)=PIMD forces on atoms in each cell
  mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
  natom=number of atoms
  trotter=Trotter number
  xcart(3,natom,trotter)=cartesian coordinates of atoms in each cell at t
  xcart_prev(3,natom,trotter)=cartesian coordinates of atoms in each cell at t-dt

OUTPUT

  xcart_next(3,natom,trotter)=cartesian coordinates of atoms in each cell at t+dt

SIDE EFFECTS

SOURCE

1820 subroutine pimd_predict_verlet(dtion,forces,mass,natom,trotter,xcart,xcart_next,xcart_prev)
1821 
1822 !Arguments ------------------------------------
1823 !scalars
1824  integer,intent(in) :: natom,trotter
1825  real(dp),intent(in) :: dtion
1826 !arrays
1827  real(dp),intent(in) :: forces(3,natom,trotter)
1828  real(dp),intent(in) :: xcart(3,natom,trotter),xcart_prev(3,natom,trotter)
1829  real(dp),intent(in) :: mass(:,:)
1830  real(dp),intent(out) :: xcart_next(3,natom,trotter)
1831 !Local variables-------------------------------
1832 !scalars
1833  integer :: iatom,ii,iimage,imass,natom_mass,nmass
1834  character(len=500) :: msg
1835 !arrays
1836 
1837 !************************************************************************
1838 
1839  natom_mass=size(mass,1);nmass=size(mass,2)
1840  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
1841    msg='Wrong dimensions for array mass !'
1842    ABI_BUG(msg)
1843  end if
1844 
1845  do iimage=1,trotter
1846    imass=min(nmass,iimage)
1847    do iatom=1,natom
1848      do ii=1,3
1849        xcart_next(ii,iatom,iimage)= &
1850 &         two*xcart(ii,iatom,iimage) &
1851 &       - xcart_prev(ii,iatom,iimage) &
1852 &       + dtion*dtion*forces(ii,iatom,iimage)/mass(iatom,imass)
1853      end do
1854    end do
1855  end do
1856 
1857 end subroutine pimd_predict_verlet

m_pimd/pimd_print [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_print

FUNCTION

  Print out results related to PIMD (for given time step)

INPUTS

  constraint=type of constraint eventually applied (on a reaction coordinate)
  constraint_output(2)=several (real) data to be output when a constraint has been applied
  eharm=harmonic energy
  eharm_virial=harmonic energy from virial
  epot=potential energy
  forces(3,natom,trotter)=forces on atoms in each cell
  inertmass(natom)=inertial masses of atoms
  irestart=1 if this is a restart
  itimimage=index of time step
  kt=Thermal energy K_b.T
  natom=number of atoms
  optcell=option for cell evolution
  prtstress=flag for stress tensor printing
  prtvolimg=printing volume
  rprim(3,3)=dimensionless real space primitive translations
  stress(3,3,3)=stress tensor (first dimension corresponds to 3 different estimators)
  temperature1,temperature2=temperatures at t and t+dt
  traj_unit=flag activating printing of the trajectory in an external file\
            if >0, indicates the unit number of the trajectory file
  trotter=Trotter number
  vel(3,natom,trotter)=velocities of atoms in in each cell
  vel_cell(3,3)= velocities of cell parameters (time derivative of rprimd)
  xcart(3,natom,trotter)=cartesian coordinates of atoms in in each cell
  xred(3,natom,trotter)=reduced coordinates of atoms in in each cell

OUTPUT

  -- only printing --

SIDE EFFECTS

SOURCE

572 subroutine pimd_print(constraint,constraint_output,eharm,eharm_virial,epot,&
573 &          forces,inertmass,irestart,itimimage,kt,natom,optcell,prtstress,&
574 &          prtvolimg,rprimd,stress,temperature1,temperature2,traj_unit,&
575 &          trotter,vel,vel_cell,xcart,xred)
576 
577 !Arguments ------------------------------------
578 !scalars
579  integer,intent(in) :: constraint,irestart,itimimage,natom,optcell
580  integer,intent(in) :: prtstress,prtvolimg,traj_unit,trotter
581  real(dp),intent(in) :: eharm,eharm_virial,epot,kt,temperature1,temperature2
582 !arrays
583  real(dp),intent(in) :: constraint_output(2)
584  real(dp),intent(in) :: forces(3,natom,trotter),inertmass(natom)
585  real(dp),intent(in) :: rprimd(3,3),stress(3,3,3),vel(3,natom,trotter),vel_cell(3,3)
586  real(dp),intent(in) :: xcart(3,natom,trotter),xred(3,natom,trotter)
587 !Local variables-------------------------------
588 !scalars
589  integer :: iatom,ii,iimage
590  real(dp) :: mtot
591  character(len=500) :: msg
592 !arrays
593  real(dp) :: acell(3),cdm_cart(3),cdm_red(3),forcetot(3),rprim(3,3)
594  real(dp),allocatable :: centroid(:,:),qudeloc(:)
595 
596 !************************************************************************
597 
598 !Temperature
599 if(itimimage==1)then
600    msg=ch10
601    if(mod(irestart,10)==0) then
602      write(msg,'(2a)') ch10,' This is a PIMD calculation from scratch'
603    else if (mod(irestart,10)==1) then
604      write(msg,'(2a)') ch10,' This is a RESTART calculation'
605    end if
606    call wrtout(ab_out,msg,'COLL')
607    call wrtout(std_out,msg,'COLL')
608    write(msg,'(a,f12.5,a)') &
609 &   ' In the initial configuration, the temperature is ',temperature1,' K'
610    call wrtout(ab_out,msg,'COLL')
611    call wrtout(std_out,msg,'COLL')
612  end if
613  write(msg,'(2a,i5,a,f12.5,a)') ch10,&
614 &  ' At PIMD time step ',itimimage,', the temperature is',temperature2,' K'
615  call wrtout(ab_out,msg,'COLL')
616  call wrtout(std_out,msg,'COLL')
617 
618 !Energies
619  write(msg,'(4a,f18.9,a,a,a,f18.9,a,a)') ch10,&
620 &  ' Energy:',ch10, &
621 &  '   Internal energy (PRIMITIVE estimator) =',onehalf*dble(natom*trotter)*kt-eharm+epot ,' Ha',ch10, &
622 &  '   Internal energy (VIRIAL    estimator) =',onehalf*dble(natom)*kt+eharm_virial+epot,' Ha',ch10
623  call wrtout(ab_out,msg,'COLL')
624  call wrtout(std_out,msg,'COLL')
625 
626 !Stress tensor and pressure
627  write(msg,'(2a,3(2a,3f18.9))') ch10,&
628 &   ' Stress tensor from PRIMITIVE estimator (Ha/Bohr^3):',ch10, &
629 &   '   ',stress(1,1,1),stress(1,1,2),stress(1,1,3),ch10, &
630 &   '   ',stress(1,2,1),stress(1,2,2),stress(1,2,3),ch10, &
631 &   '   ',stress(1,3,1),stress(1,3,2),stress(1,3,3)
632  if (prtstress==1) then
633    call wrtout(ab_out,msg,'COLL')
634  end if
635  call wrtout(std_out,msg,'COLL')
636  write(msg,'(a,f18.9,a)') ' Pressure (primitive estimator) =', &
637 &  -third*(stress(1,1,1)+stress(1,2,2)+stress(1,3,3))*HaBohr3_GPa,' GPa'
638  call wrtout(ab_out,msg,'COLL')
639  call wrtout(std_out,msg,'COLL')
640 
641 !Data related to constraint eventually applied
642  if (constraint/=0) then
643    if (constraint==1) write(msg,'(2a)') ch10,' Blue Moon Ensemble method is activated:'
644    call wrtout(ab_out,msg,'COLL')
645    call wrtout(std_out,msg,'COLL')
646    write(msg,'(a,f18.10,2a,f18.10)') &
647 &     '  - Reaction coordinate =',constraint_output(1),ch10,&
648 &     '  - Instantaneous force on the reaction coord. =',constraint_output(2)
649    call wrtout(ab_out,msg,'COLL')
650    call wrtout(std_out,msg,'COLL')
651  end if
652 
653 !Total force
654  if (prtvolimg<=1) then
655    forcetot=zero
656    do iimage=1,trotter
657      do iatom=1,natom
658        do ii=1,3
659          forcetot(ii)=forcetot(ii)+forces(ii,iatom,iimage)
660        end do
661      end do
662    end do
663    write(msg,'(2a,3f18.10)') ch10,' Total force=',forcetot(1:3)
664    call wrtout(std_out,msg,'COLL')
665  end if
666 
667 !position of mass center
668  if (prtvolimg<=1) then
669    mtot=zero;cdm_cart=zero;cdm_red=zero
670    do iimage=1,trotter
671      do iatom=1,natom
672        cdm_cart(:)=cdm_cart(:)+inertmass(iatom)*xcart(:,iatom,iimage)
673        cdm_red (:)=cdm_red (:)+inertmass(iatom)*xred (:,iatom,iimage)
674        mtot=mtot+inertmass(iatom)
675      end do
676    end do
677    cdm_cart=cdm_cart/mtot
678    cdm_red =cdm_red /mtot
679    write(msg,'(3a,3x,3f18.10,3a,3x,3f18.10)') ch10,&
680 &    ' Center of mass, in cartes. coordinates :',ch10,cdm_cart(:),ch10,&
681 &    ' Center of mass, in reduced coordinates :',ch10,cdm_red(:)
682    call wrtout(std_out,msg,'COLL')
683    call wrtout(ab_out,msg,'COLL')
684  end if
685 
686 !Positions
687  write(msg,'(2a)') ch10,' Atomic positions:'
688  call wrtout(std_out,msg,'COLL')
689  call wrtout(ab_out,msg,'COLL')
690  do iimage=1,trotter
691    select case(iimage)
692    case(1)
693     write(msg,'(a)') ' xred'
694    case(2,3,4,5,6,7,8,9)
695      write(msg,'(a,i1,a)') ' xred_',iimage,'img'
696    case(10:99)
697      write(msg,'(a,i2,a)') ' xred_',iimage,'img'
698    case default
699      write(msg,'(a,i3,a)') ' xred_',iimage,'img'
700    end select
701    call wrtout(std_out,msg,'COLL')
702    call wrtout(ab_out,msg,'COLL')
703    if (traj_unit>0) then
704      call wrtout(traj_unit,msg,'COLL')
705    end if
706    do iatom=1,natom
707      write(msg,'(3f18.10)') xred(1:3,iatom,iimage)
708      call wrtout(std_out,msg,'COLL')
709      call wrtout(ab_out,msg,'COLL')
710      if (traj_unit>0) then
711        call wrtout(traj_unit,msg,'COLL')
712      end if
713    end do
714  end do
715 
716 !Velocities
717  write(msg,'(2a)') ch10,' Velocities:'
718  call wrtout(std_out,msg,'COLL')
719  call wrtout(ab_out,msg,'COLL')
720  do iimage=1,trotter
721    select case(iimage)
722    case(1)
723      write(msg,'(a)') ' vel'
724    case(2,3,4,5,6,7,8,9)
725      write(msg,'(a,i1,a)') ' vel_',iimage,'img'
726    case(10:99)
727      write(msg,'(a,i2,a)') ' vel_',iimage,'img'
728    case default
729      write(msg,'(a,i3,a)') ' vel_',iimage,'img'
730    end select
731    call wrtout(std_out,msg,'COLL')
732    call wrtout(ab_out,msg,'COLL')
733    if (traj_unit>0) then
734      call wrtout(traj_unit,msg,'COLL')
735    end if
736    do iatom=1,natom
737      write(msg,'(3f18.10)') vel(1:3,iatom,iimage)
738      call wrtout(std_out,msg,'COLL')
739      call wrtout(ab_out,msg,'COLL')
740      if (traj_unit>0) then
741        call wrtout(traj_unit,msg,'COLL')
742      end if
743    end do
744  end do
745 
746  if (optcell>0) then
747 
748    call mkradim(acell,rprim,rprimd)
749 
750 !  Time derivative of rprimd
751    write(msg,'(2a)') ch10,' Time derivative of rprimd:'
752    call wrtout(std_out,msg,'COLL')
753    call wrtout(ab_out,msg,'COLL')
754    write(msg,'(2a,3(3f18.10))') ' vel_cell',ch10,vel_cell(:,:)
755    call wrtout(std_out,msg,'COLL')
756    call wrtout(ab_out,msg,'COLL')
757    if (traj_unit>0) then
758      call wrtout(traj_unit,msg,'COLL')
759    end if
760 
761 !  rprimd
762    write(msg,'(2a)') ch10,' Cell parameters:'
763    call wrtout(std_out,msg,'COLL')
764    call wrtout(ab_out,msg,'COLL')
765    write(msg,'(2a,3(3f18.10),3a,3f18.10)') ' rprim',ch10,rprim(:,:),ch10,&
766 &                                          ' acell',ch10,acell(:)
767    call wrtout(std_out,msg,'COLL')
768    call wrtout(ab_out,msg,'COLL')
769    if (traj_unit>0) then
770      call wrtout(traj_unit,msg,'COLL')
771    end if
772 
773  end if
774 
775 !Centroids and wave-packet spatial spreads
776  ABI_MALLOC(centroid,(3,natom))
777  ABI_MALLOC(qudeloc,(natom))
778  centroid=zero;qudeloc=zero
779  do iimage=1,trotter
780    do iatom=1,natom
781      do ii=1,3
782        centroid(ii,iatom)=centroid(ii,iatom)+xcart(ii,iatom,iimage)
783      end do
784    end do
785  end do
786  centroid=centroid/dble(trotter)
787  do iimage=1,trotter
788    do iatom=1,natom
789      do ii=1,3
790        qudeloc(iatom)=qudeloc(iatom)+((xcart(ii,iatom,iimage)-centroid(ii,iatom))**2)
791      end do
792    end do
793  end do
794  qudeloc(:)=sqrt(qudeloc(:)/dble(trotter))
795  write(msg,'(4a)') ch10,' Centroids and wave-packet spatial spreads (cart. coord.):',ch10,&
796 &  ' iat        centroid_x        centroid_y        centroid_z    spatial_spread'
797  call wrtout(std_out,msg,'COLL')
798  do iatom=1,natom
799    write(msg,'(i4,4f18.10)') iatom,centroid(1:3,iatom),qudeloc(iatom)
800    call wrtout(std_out,msg,'COLL')
801  end do
802  ABI_FREE(centroid)
803  ABI_FREE(qudeloc)
804 
805 !Fake statement
806  return;ii=prtvolimg
807 
808 end subroutine pimd_print

m_pimd/pimd_skip_qtb [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_skip_qtb

FUNCTION

  Only relevant in case of PI-QTB:
  Skip a line in a QTB random force file

INPUTS

  pimd_param=datastructure of type pimd_type.
             several parameters for Path-Integral MD.

OUTPUT

SOURCE

398 subroutine pimd_skip_qtb(pimd_param)
399 
400 !Arguments ------------------------------------
401 !scalars
402  type(pimd_type),intent(in) :: pimd_param
403 !arrays
404 !Local variables-------------------------------
405 !scalars
406  character(len=500) :: msg
407 !arrays
408 
409 !************************************************************************
410 
411  if (pimd_param%use_qtb==0) return
412 
413  if (pimd_param%qtb_file_unit<0) then
414    msg='QTB forces file unit should be positive!'
415    ABI_BUG(msg)
416  end if
417 
418 !Skip one line QTB random forces file
419  read(pimd_param%qtb_file_unit)
420 
421 end subroutine pimd_skip_qtb

m_pimd/pimd_stresses [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_stresses

FUNCTION

  In the case od PIMD, compute the pressure tensor from virial theorem

INPUTS

  mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
  natom=number of atoms
  quantummass(natom)=quantum masses of atoms
  stressin(3,3,trotter)=electronic stress tensor for each image
  temperature=temperature (could be instantaneous temp. or thermostat temp.)
  temperature_therm=thermostat temperature
  trotter=Trotter number
  vel(3,natom,trotter)=velocities of atoms in each cell
  volume=volume of each cell (common to all cells)
  xcart(3,natom,trotter)=cartesian coordinates of atoms in each cell at t
  temperature=thermostat temperature

OUTPUT

  stress_pimd(3,3,3)=stress tensor for PIMD
                     First dimension (3) corresponds to 3 different pressure estimators

SIDE EFFECTS

SOURCE

1565 subroutine pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,&
1566 &                        temperature,temperature_therm,trotter,vel,volume,xcart)
1567 
1568 !Arguments ------------------------------------
1569 !scalars
1570  integer,intent(in) :: natom,trotter
1571  real(dp),intent(in) :: temperature,temperature_therm,volume
1572 !arrays
1573  real(dp),intent(in) :: vel(3,natom,trotter),xcart(3,natom,trotter)
1574  real(dp),intent(in) :: mass(:,:),stressin(3,3,trotter),quantummass(natom)
1575  real(dp),intent(out) :: stress_pimd(3,3,3)
1576 !Local variables-------------------------------
1577 !scalars
1578  integer :: iatom,ii,iimage,imass,jj,natom_mass,nmass,iimagep
1579  real(dp) :: stress_tmp(3,3,3),omega2,kt,kt_therm
1580  character(len=500) :: msg
1581 
1582 !************************************************************************
1583 
1584  natom_mass=size(mass,1);nmass=size(mass,2)
1585  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
1586    msg='Wrong dimensions for array mass !'
1587    ABI_BUG(msg)
1588  end if
1589 
1590  kt=temperature*kb_HaK
1591  kt_therm=temperature_therm*kb_HaK
1592  stress_pimd=zero
1593 
1594 !I-PRIMITIVE ESTIMATOR
1595 !1-Kinetic part
1596  do ii=1,3
1597    stress_pimd(1,ii,ii)=dble(natom)*dble(trotter)*kt/volume
1598  end do
1599 
1600 !2-Potential part
1601  do iimage=1,trotter
1602    do ii=1,3
1603      do jj=1,3
1604        stress_pimd(1,ii,jj)=stress_pimd(1,ii,jj)-stressin(ii,jj,iimage)/dble(trotter)
1605        !minus to convert stress into pressure
1606      end do
1607    end do
1608  end do
1609 
1610 !3-Contribution from springs
1611  omega2=dble(trotter)*kt_therm*kt_therm
1612  do iimage=1,trotter
1613    iimagep=iimage+1
1614    if(iimage==trotter) then
1615      iimagep=1
1616    end if
1617    do iatom=1,natom
1618      do ii=1,3
1619        do jj=1,3
1620          stress_pimd(1,ii,jj)=stress_pimd(1,ii,jj)-quantummass(iatom)*omega2* &
1621 &        (xcart(ii,iatom,iimagep)-xcart(ii,iatom,iimage))*  &
1622 &        (xcart(jj,iatom,iimagep)-xcart(jj,iatom,iimage))/volume
1623        end do
1624      end do
1625    end do
1626  end do
1627 
1628 !II-Average of classical pressures
1629 !1-Kinetic part
1630  do iimage=1,trotter
1631   imass=min(nmass,iimage)
1632   do iatom=1,natom
1633     do ii=1,3
1634       do jj=1,3
1635         stress_pimd(2,ii,jj)=stress_pimd(2,ii,jj)+ &
1636 &          mass(iatom,imass)*vel(ii,iatom,iimage)*vel(jj,iatom,iimage)/volume
1637       end do
1638     end do
1639   end do
1640  end do
1641 
1642 !2-Contribution from electronic stress
1643  do iimage=1,trotter
1644    do ii=1,3
1645      do jj=1,3
1646        stress_pimd(2,ii,jj)=stress_pimd(2,ii,jj)-stressin(ii,jj,iimage)
1647      end do
1648    end do
1649  end do
1650  do ii=1,3
1651    do jj=1,3
1652      stress_pimd(2,ii,jj)=stress_pimd(2,ii,jj)/dble(trotter)
1653    end do
1654  end do
1655 
1656 !III-pressure from VIRIAL estimator: stress_pimd(3,:,:)
1657 !1-kinetic part
1658 ! do ii=1,3
1659 !   stress_pimd(3,ii,ii)=dfloat(natom)*kt/volume
1660 ! end do
1661 
1662 !Symmetrize internal pressure
1663  stress_tmp=stress_pimd
1664  stress_pimd(:,2,1)=half*(stress_tmp(:,1,2)+stress_tmp(:,2,1))
1665  stress_pimd(:,1,2)=half*(stress_tmp(:,1,2)+stress_tmp(:,2,1))
1666  stress_pimd(:,3,1)=half*(stress_tmp(:,1,3)+stress_tmp(:,3,1))
1667  stress_pimd(:,1,3)=half*(stress_tmp(:,1,3)+stress_tmp(:,3,1))
1668  stress_pimd(:,2,3)=half*(stress_tmp(:,3,2)+stress_tmp(:,2,3))
1669  stress_pimd(:,3,2)=half*(stress_tmp(:,3,2)+stress_tmp(:,2,3))
1670 
1671 end subroutine pimd_stresses

m_pimd/pimd_temperature [ Functions ]

[ Top ] [ m_pimd ] [ Functions ]

NAME

  pimd_temperature

FUNCTION

  Compute temperature from velocities and masses

INPUTS

  mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter)
  vel(3,natom,nimage)=velocities for each image of the cell

OUTPUT

  pimd_temperature=temperature (from all images of the cell)

SIDE EFFECTS

SOURCE

489 function pimd_temperature(mass,vel)
490 
491 !Arguments ------------------------------------
492 !scalars
493  real(dp) :: pimd_temperature
494 !arrays
495  real(dp),intent(in) :: mass(:,:),vel(:,:,:)
496 !Local variables-------------------------------
497 !scalars
498  integer :: iatom,idir,iimage,imass,natom,natom_mass,ndir,nimage,nmass
499  real(dp) :: v2
500  character(len=500) :: msg
501 !arrays
502 
503 !************************************************************************
504 
505  ndir=size(vel,1);natom=size(vel,2);nimage=size(vel,3)
506  natom_mass=size(mass,1);nmass=size(mass,2)
507  if (ndir/=3.or.natom<=0.or.nimage<=0) then
508    msg='Wrong sizes for vel array !'
509    ABI_BUG(msg)
510  end if
511  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=nimage)) then
512    msg='Wrong dimensions for array mass !'
513    ABI_BUG(msg)
514  end if
515 
516  v2=zero
517  do iimage=1,nimage
518    imass=min(nmass,iimage)
519    do iatom=1,natom
520      do idir=1,3
521        v2=v2+vel(idir,iatom,iimage)*vel(idir,iatom,iimage)*mass(iatom,imass)
522      end do
523    end do
524  end do
525  pimd_temperature=v2/(dble(3*natom*nimage)*kb_HaK)
526 
527 end function pimd_temperature

m_pimd/pimd_type [ Types ]

[ Top ] [ m_pimd ] [ Types ]

NAME

 pimd_type

FUNCTION

 Datatype with the variables required to perform PIMD

NOTES

SOURCE

 80  type,public :: pimd_type
 81 ! Scalars
 82   integer  :: adpimd
 83   integer  :: constraint
 84   integer  :: irandom
 85   integer  :: nnos
 86   integer  :: ntypat
 87   integer  :: optcell
 88   integer  :: pitransform
 89   integer  :: traj_unit
 90   integer  :: use_qtb
 91   integer  :: qtb_file_unit
 92   real(dp) :: adpimd_gamma
 93   real(dp) :: vis
 94   real(dp) :: bmass
 95   real(dp) :: dtion
 96   real(dp) :: friction
 97 ! Arrays
 98   integer ,pointer  :: typat(:)      ! This pointer is associated with dtset%typat
 99   real(dp),pointer :: amu(:)         ! This pointer is associated with dtset%%amu_orig(:,1)
100   real(dp),pointer :: mdtemp(:)      ! This pointer is associated with dtset%mdtemp
101   real(dp),pointer :: pimass(:)      ! This pointer is associated with dtset%pimass
102   real(dp),pointer :: qmass(:)       ! This pointer is associated with dtset%qmass
103   real(dp),pointer :: strtarget(:)   ! This pointer is associated with dtset%strtarget
104   real(dp),pointer :: wtatcon(:,:,:) ! This pointer is associated with dtset%wtatcon
105   real(dp),allocatable :: zeta_prev(:,:,:,:)
106   real(dp),allocatable :: zeta     (:,:,:,:)
107   real(dp),allocatable :: zeta_next(:,:,:,:)
108   real(dp),allocatable :: dzeta    (:,:,:,:)
109  end type pimd_type