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-2018 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 .

PARENTS

CHILDREN

SOURCE

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

PARENTS

      pimd_langevin_npt,pimd_langevin_nvt,pimd_nosehoover_nvt

CHILDREN

SOURCE

2709 subroutine pimd_apply_constraint(constraint,constraint_output,forces,mass,natom,&
2710 &                                trotter,wtatcon,xcart)
2711 
2712 
2713 !This section has been created automatically by the script Abilint (TD).
2714 !Do not modify the following lines by hand.
2715 #undef ABI_FUNC
2716 #define ABI_FUNC 'pimd_apply_constraint'
2717 !End of the abilint section
2718 
2719  implicit none
2720 
2721 !Arguments ------------------------------------
2722 !scalars
2723  integer,intent(in) :: constraint,natom,trotter
2724 !arrays
2725  real(dp),intent(in) :: mass(natom),wtatcon(3,natom),xcart(3,natom,trotter)
2726  real(dp),intent(out) :: constraint_output(2)
2727  real(dp),intent(inout) :: forces(3,natom,trotter)
2728 !Local variables-------------------------------
2729 !scalars
2730  integer :: iatom,ii,iimage
2731  real(dp) :: af,lambda_cst,masstot,one_over_trotter,xcart_centroid,zz
2732  !character(len=500) :: msg
2733 !arrays
2734  real(dp) :: force_centroid(3),lambda_com(3),mat(3,3),matinv(3,3),vec(3),weightsum(3)
2735 
2736 !************************************************************************
2737 
2738 
2739  select case(constraint)
2740 
2741 !=== No constraint =======================================================
2742  case(0)
2743 
2744    constraint_output(:)=zero
2745    return
2746 
2747 !=== Linear combination of centroid coordinates ==========================
2748  case(1)
2749 
2750 !  Some useful quantities
2751    masstot=sum(mass)*dble(trotter)
2752    weightsum(:)=sum(wtatcon,dim=2)*dble(trotter)
2753    zz=zero;af=zero;force_centroid=zero
2754    do iimage=1,trotter
2755      do iatom=1,natom
2756        do ii=1,3
2757          force_centroid(ii)=force_centroid(ii)+forces(ii,iatom,iimage)
2758          af=af+wtatcon(ii,iatom)*forces(ii,iatom,iimage)/mass(iatom)
2759          zz=zz+wtatcon(ii,iatom)**2/mass(iatom)
2760        end do
2761      end do
2762    end do
2763    vec(:)=force_centroid(:)-weightsum(:)*af/zz
2764    do ii=1,3
2765      mat(:,ii)=-weightsum(:)*weightsum(ii)/zz
2766      mat(ii,ii)=mat(ii,ii)+masstot
2767    end do
2768    call matr3inv(mat,matinv)
2769 
2770    !Calculation of a Lagrange multipliers:
2771    ! lambda_cst: to apply the constraint
2772    ! lambda_com: to maintain the position of the center of mass
2773    lambda_com(:)=matmul(matinv,vec)
2774    lambda_cst=(af-dot_product(weightsum,lambda_com))*dble(trotter)/zz
2775 
2776    !Modification of forces
2777    one_over_trotter=one/dble(trotter)
2778    do iimage=1,trotter
2779      do iatom=1,natom
2780        do ii=1,3
2781          forces(ii,iatom,iimage)=forces(ii,iatom,iimage) &
2782 &                               -lambda_cst*wtatcon(ii,iatom)*one_over_trotter &
2783 &                               -lambda_com(ii)*mass(iatom)
2784        end do
2785      end do
2786    end do
2787 
2788    !Computation of relevant outputs
2789    constraint_output(:)=zero
2790    !1-Reaction coordinate
2791    do iatom=1,natom
2792      do ii=1,3
2793        xcart_centroid=zero
2794        do iimage=1,trotter
2795          xcart_centroid=xcart_centroid+xcart(ii,iatom,iimage)
2796        end do
2797        constraint_output(1)=constraint_output(1)+xcart_centroid*wtatcon(ii,iatom)
2798      end do
2799    end do
2800    constraint_output(1)=constraint_output(1)/dble(trotter)
2801    !2-Force on reaction coordinate
2802    constraint_output(2)=-lambda_cst
2803 
2804  end select
2805 
2806 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

PARENTS

      pimd_langevin_nvt,pimd_nosehoover_nvt

CHILDREN

SOURCE

2378 subroutine pimd_coord_transform(array,ioption,natom,transform,trotter)
2379 
2380 
2381 !This section has been created automatically by the script Abilint (TD).
2382 !Do not modify the following lines by hand.
2383 #undef ABI_FUNC
2384 #define ABI_FUNC 'pimd_coord_transform'
2385 !End of the abilint section
2386 
2387  implicit none
2388 
2389 !Arguments ------------------------------------
2390 !scalars
2391  integer,intent(in) :: ioption,natom,transform,trotter
2392 !arrays
2393  real(dp),intent(inout) :: array(3,natom,trotter)
2394 !Local variables-------------------------------
2395 !scalars
2396  integer :: iatom,ii,iimage,iimagem,iimagep,iimage2,ll
2397 !arrays
2398  real(dp),allocatable :: array_temp(:,:,:)
2399  real(dp),allocatable :: nrm(:,:)
2400 
2401 !************************************************************************
2402 
2403 !=== No transformation ===================================================
2404  if (transform==0) then
2405 
2406    return
2407 
2408 !=== Normal mode transformation ==========================================
2409  else if (transform==1) then
2410 
2411 !  ---------------- From primitive to transformed coordinates ------------
2412    if (ioption==+1) then
2413      ABI_ALLOCATE(array_temp,(3,natom,trotter))
2414 
2415      array_temp(:,:,1)=zero
2416      do iimage=1,trotter
2417        array_temp(:,:,1)=array_temp(:,:,1)+array(:,:,iimage)
2418      end do
2419 
2420      do ll=2,trotter/2
2421        array_temp(:,:,2*ll-2)=zero
2422        array_temp(:,:,2*ll-1)=zero
2423        do iimage=1,trotter
2424          array_temp(:,:,2*ll-2)=array_temp(:,:,2*ll-2)+array(:,:,iimage)* &
2425 &                              cos(two_pi*dble((iimage-1)*(ll-1))/dble(trotter))
2426          array_temp(:,:,2*ll-1)=array_temp(:,:,2*ll-1)-array(:,:,iimage)* &
2427 &                              sin(two_pi*dble((iimage-1)*(ll-1))/dble(trotter))
2428        end do
2429      end do
2430 
2431      array_temp(:,:,trotter)=zero
2432      do iimage=1,trotter
2433        array_temp(:,:,trotter)=array_temp(:,:,trotter)+&
2434 &                              array(:,:,iimage)*dble((-1)**(iimage-1))
2435      end do
2436 
2437      array=array_temp/dble(trotter)
2438 
2439      ABI_DEALLOCATE(array_temp)
2440 
2441 !  ---------------- From transformed to primitive coordinates ------------
2442    else if (ioption==-1) then
2443 
2444     ABI_ALLOCATE(array_temp,(3,natom,trotter))  !real part
2445     ABI_ALLOCATE(nrm,(trotter,trotter))  !real part
2446 
2447    do iimage=1,trotter
2448      nrm(iimage,1)=one
2449      nrm(iimage,trotter)=dble((-1)**(iimage-1))
2450    end do
2451 
2452    do ll=2,trotter/2
2453      do iimage=1,trotter
2454        nrm(iimage,2*ll-2)= two*cos(two_pi*dble((ll-1)* &
2455 &                        (iimage-1))/dble(trotter))
2456        nrm(iimage,2*ll-1)=-two*sin(two_pi*dble((ll-1)* &
2457 &                        (iimage-1))/dble(trotter))
2458      end do
2459    end do
2460 
2461     do iimage=1,trotter
2462       array_temp(:,:,iimage)=zero
2463       do ll=1,trotter
2464         array_temp(:,:,iimage)=array_temp(:,:,iimage)+nrm(iimage,ll)*array(:,:,ll)
2465       end do
2466     end do
2467     array=array_temp
2468 
2469     ABI_DEALLOCATE(array_temp)
2470     ABI_DEALLOCATE(nrm)
2471 
2472    end if ! ioption
2473 
2474 !=== Staging transformation ==============================================
2475  else if (transform==2) then
2476 
2477 !  ---------------- From primitive to transformed coordinates ------------
2478    if (ioption==+1) then
2479 
2480      ABI_ALLOCATE(array_temp,(3,natom,trotter))
2481      array_temp=zero
2482      do iimage=1,trotter
2483        iimagep=iimage+1;if(iimage==trotter) iimagep=1
2484        iimagem=iimage-1;if(iimage==1) iimagem=trotter
2485        do iatom=1,natom
2486          do ii=1,3
2487            array_temp(ii,iatom,iimage)=(dble(iimagem)*array(ii,iatom,iimagep)+array(ii,iatom,1))/dble(iimage)
2488          end do
2489        end do
2490      end do
2491      if (trotter>1) then
2492        do iimage=2,trotter
2493          do iatom=1,natom
2494            do ii=1,3
2495              array(ii,iatom,iimage)=array(ii,iatom,iimage)-array_temp(ii,iatom,iimage)
2496            end do
2497          end do
2498        end do
2499      end if
2500      ABI_DEALLOCATE(array_temp)
2501 
2502 !  ---------------- From transformed to primitive coordinates ------------
2503    else if (ioption==-1) then
2504 
2505      ABI_ALLOCATE(array_temp,(3,natom,trotter))
2506      array_temp=zero
2507      do iimage=1,trotter
2508        do iatom=1,natom
2509          do ii=1,3
2510            array_temp(ii,iatom,iimage)=array(ii,iatom,1)
2511          end do
2512        end do
2513      end do
2514      if (trotter>1) then
2515        do iimage=2,trotter
2516          do iatom=1,natom
2517            do ii=1,3
2518              do iimage2=iimage,trotter
2519                array_temp(ii,iatom,iimage)=array_temp(ii,iatom,iimage) &
2520 &                     +array(ii,iatom,iimage2)*(dble(iimage-1))/(dble(iimage2-1))
2521              end do
2522            end do
2523          end do
2524        end do
2525      end if
2526      array(:,:,:)=array_temp(:,:,:)
2527      ABI_DEALLOCATE(array_temp)
2528 
2529    end if ! ioption
2530 
2531  end if ! transform
2532 
2533 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.

PARENTS

      gstateimg

CHILDREN

SOURCE

309 subroutine pimd_destroy(pimd_param)
310 
311 
312 !This section has been created automatically by the script Abilint (TD).
313 !Do not modify the following lines by hand.
314 #undef ABI_FUNC
315 #define ABI_FUNC 'pimd_destroy'
316 !End of the abilint section
317 
318  implicit none
319 
320 !Arguments ------------------------------------
321 !scalars
322  type(pimd_type),intent(inout) :: pimd_param
323 !arrays
324 !Local variables-------------------------------
325 !scalars
326  integer :: ierr
327  character(len=100) :: msg
328 !arrays
329 
330 !************************************************************************
331 
332  if (allocated(pimd_param%zeta_prev)) then
333    ABI_DEALLOCATE(pimd_param%zeta_prev)
334  end if
335  if (allocated(pimd_param%zeta)) then
336    ABI_DEALLOCATE(pimd_param%zeta)
337  end if
338  if (allocated(pimd_param%zeta_next)) then
339    ABI_DEALLOCATE(pimd_param%zeta_next)
340  end if
341  if (allocated(pimd_param%dzeta)) then
342    ABI_DEALLOCATE(pimd_param%dzeta)
343  end if
344 
345  if (pimd_param%qtb_file_unit>0) then
346    if (is_open(pimd_param%qtb_file_unit)) then
347      ierr=close_unit(pimd_param%qtb_file_unit,msg)
348    end if
349  end if
350 
351  if (pimd_param%traj_unit>0) then
352    if (is_open(pimd_param%traj_unit)) then
353      ierr=close_unit(pimd_param%traj_unit,msg)
354    end if
355  end if
356 
357  call pimd_nullify(pimd_param)
358 
359 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

PARENTS

CHILDREN

SOURCE

1954 function pimd_diff_stress(stress_pimd,stress_target)
1955 
1956 
1957 !This section has been created automatically by the script Abilint (TD).
1958 !Do not modify the following lines by hand.
1959 #undef ABI_FUNC
1960 #define ABI_FUNC 'pimd_diff_stress'
1961 !End of the abilint section
1962 
1963  implicit none
1964 
1965 !Arguments ------------------------------------
1966 !scalars
1967 !arrays
1968  real(dp),intent(in) :: stress_pimd(3,3,3),stress_target(6)
1969  real(dp) :: pimd_diff_stress(3,3)
1970 !Local variables-------------------------------
1971 !scalars
1972  integer :: ii,jj
1973 !arrays
1974  real(dp) :: stress_pimd2(3,3)
1975 
1976 !************************************************************************
1977 
1978 !Choice: the primitive estimator for pressure is chosen
1979  do ii=1,3
1980    do jj=1,3
1981      stress_pimd2(ii,jj)=stress_pimd(1,ii,jj)
1982    end do
1983  end do
1984 
1985 !+stress_target instead of - because it is translated from stress to pressure tensor
1986  pimd_diff_stress(1,1)=stress_pimd2(1,1)+stress_target(1)
1987  pimd_diff_stress(2,2)=stress_pimd2(2,2)+stress_target(2)
1988  pimd_diff_stress(3,3)=stress_pimd2(3,3)+stress_target(3)
1989  pimd_diff_stress(2,3)=stress_pimd2(2,3)+stress_target(4)
1990  pimd_diff_stress(3,2)=stress_pimd2(3,2)+stress_target(4)
1991  pimd_diff_stress(1,3)=stress_pimd2(1,3)+stress_target(5)
1992  pimd_diff_stress(3,1)=stress_pimd2(3,1)+stress_target(5)
1993  pimd_diff_stress(1,2)=stress_pimd2(1,2)+stress_target(6)
1994  pimd_diff_stress(2,1)=stress_pimd2(2,1)+stress_target(6)
1995 
1996 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

PARENTS

      pimd_langevin_npt,pimd_langevin_nvt,pimd_nosehoover_npt
      pimd_nosehoover_nvt

CHILDREN

SOURCE

1448 subroutine pimd_energies(eharm,eharm_virial,epot,etotal_img,forces,natom,spring,trotter,xcart)
1449 
1450 
1451 !This section has been created automatically by the script Abilint (TD).
1452 !Do not modify the following lines by hand.
1453 #undef ABI_FUNC
1454 #define ABI_FUNC 'pimd_energies'
1455 !End of the abilint section
1456 
1457  implicit none
1458 
1459 !Arguments ------------------------------------
1460 !scalars
1461  integer,intent(in) :: natom,trotter
1462  real(dp),intent(out) :: eharm,eharm_virial,epot
1463 !arrays
1464  real(dp),intent(in) :: etotal_img(trotter),forces(3,natom,trotter)
1465  real(dp),intent(in) :: xcart(3,natom,trotter)
1466  real(dp),intent(in) :: spring(natom)
1467 !Local variables-------------------------------
1468 !scalars
1469  integer :: iatom,ii,iimage,iimagep
1470 !arrays
1471  real(dp),allocatable :: centroid(:,:)
1472 
1473 !************************************************************************
1474 
1475 !Compute the centroid
1476  ABI_ALLOCATE(centroid,(3,natom))
1477  centroid=zero
1478  do iimage=1,trotter
1479    do iatom=1,natom
1480      do ii=1,3
1481        centroid(ii,iatom)=centroid(ii,iatom)+xcart(ii,iatom,iimage)
1482      end do
1483    end do
1484  end do
1485  centroid=centroid/dble(trotter)
1486 
1487 !Potential energy
1488  epot=sum(etotal_img(1:trotter))/dble(trotter)
1489 
1490 !Harmonic energy
1491  eharm=zero
1492  do iimage=1,trotter
1493    iimagep=iimage+1;if(iimage==trotter)iimagep=1
1494    do iatom=1,natom
1495      do ii=1,3
1496        eharm=eharm+half*spring(iatom)*((xcart(ii,iatom,iimagep)-xcart(ii,iatom,iimage))**2)
1497      end do
1498    end do
1499  end do
1500 
1501 !Harmonic energy from virial estimator
1502  eharm_virial=zero
1503  do iimage=1,trotter
1504    do iatom=1,natom
1505      do ii=1,3
1506        eharm_virial=eharm_virial-(xcart(ii,iatom,iimage)-centroid(ii,iatom)) &
1507 &              *forces(ii,iatom,iimage)
1508      end do
1509    end do
1510  end do
1511  eharm_virial=eharm_virial/dble(two*trotter)
1512 
1513  ABI_DEALLOCATE(centroid)
1514 
1515 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 !

PARENTS

      pimd_langevin_nvt,pimd_nosehoover_nvt

CHILDREN

SOURCE

2575 subroutine pimd_force_transform(forces,ioption,natom,transform,trotter)
2576 
2577 
2578 !This section has been created automatically by the script Abilint (TD).
2579 !Do not modify the following lines by hand.
2580 #undef ABI_FUNC
2581 #define ABI_FUNC 'pimd_force_transform'
2582 !End of the abilint section
2583 
2584  implicit none
2585 
2586 !Arguments ------------------------------------
2587 !scalars
2588  integer,intent(in) :: ioption,natom,transform,trotter
2589 !arrays
2590  real(dp),intent(inout) :: forces(3,natom,trotter)
2591 !Local variables-------------------------------
2592 !scalars
2593  integer :: iatom,ii,iimage,ll
2594  character(len=500) :: msg
2595 !arrays
2596  real(dp),allocatable :: forces_temp(:,:,:),nrm(:,:)
2597 
2598 !************************************************************************
2599 
2600  if (ioption==-1) then
2601    msg='Back transformation not implemented !'
2602    MSG_BUG(msg)
2603  end if
2604 
2605 !=== No transformation ===================================================
2606  if (transform==0) then
2607 
2608    return
2609 
2610 !=== Normal mode transformation ==========================================
2611  else if (transform==1) then
2612 
2613    !normal mode forces
2614    ABI_ALLOCATE(forces_temp,(3,natom,trotter))
2615    ABI_ALLOCATE(nrm,(trotter,trotter))
2616 
2617    do iimage=1,trotter
2618      nrm(iimage,1)=one
2619      nrm(iimage,trotter)=dble((-1)**(iimage-1))
2620    end do
2621 
2622    do ll=2,trotter/2
2623      do iimage=1,trotter
2624        nrm(iimage,2*ll-2)= two*cos(two_pi*dble((ll-1)* &
2625 &                         (iimage-1))/dble(trotter))
2626        nrm(iimage,2*ll-1)=-two*sin(two_pi*dble((ll-1)* &
2627 &                         (iimage-1))/dble(trotter))
2628      end do
2629    end do
2630 
2631    do ll=1,trotter
2632      forces_temp(:,:,ll)=zero
2633      do iimage=1,trotter
2634        forces_temp(:,:,ll)=forces_temp(:,:,ll)+nrm(iimage,ll)*forces(:,:,iimage)
2635      end do
2636    end do
2637 
2638    forces=forces_temp
2639 
2640    ABI_DEALLOCATE(forces_temp)
2641    ABI_DEALLOCATE(nrm)
2642 
2643 !=== Staging transformation ==============================================
2644  else if (transform==2) then
2645 
2646    !staging forces
2647    ABI_ALLOCATE(forces_temp,(3,natom,trotter))
2648    forces_temp=zero
2649    do iimage=1,trotter
2650      do iatom=1,natom
2651        do ii=1,3
2652          forces_temp(ii,iatom,1)=forces_temp(ii,iatom,1)+forces(ii,iatom,iimage)
2653        end do
2654      end do
2655    end do
2656    if (trotter>1) then
2657      do iimage=2,trotter
2658        do iatom=1,natom
2659          do ii=1,3
2660            forces_temp(ii,iatom,iimage)=forces(ii,iatom,iimage) &
2661 &               +forces_temp(ii,iatom,iimage-1)*(dble(iimage-2)/dble(iimage-1))
2662          end do
2663        end do
2664      end do
2665    end if
2666    forces=forces_temp
2667    ABI_DEALLOCATE(forces_temp)
2668 
2669  end if ! transform
2670 
2671 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

PARENTS

      pimd_langevin_npt,pimd_langevin_nvt,pimd_nosehoover_nvt

CHILDREN

SOURCE

1551 subroutine pimd_forces(forces,natom,spring,transform,trotter,xcart)
1552 
1553 
1554 !This section has been created automatically by the script Abilint (TD).
1555 !Do not modify the following lines by hand.
1556 #undef ABI_FUNC
1557 #define ABI_FUNC 'pimd_forces'
1558 !End of the abilint section
1559 
1560  implicit none
1561 
1562 !Arguments ------------------------------------
1563 !scalars
1564  integer,intent(in) :: natom,transform,trotter
1565 !arrays
1566  real(dp),intent(in) :: xcart(3,natom,trotter)
1567  real(dp),intent(in) :: spring(:,:)
1568  real(dp),intent(inout) :: forces(3,natom,trotter)
1569 !Local variables-------------------------------
1570 !scalars
1571  integer :: iatom,ii,iimage,iimagem,iimagep,ispring,natom_spring,nspring
1572  character(len=500) :: msg
1573 !arrays
1574 
1575 !************************************************************************
1576 
1577  natom_spring=size(spring,1);nspring=size(spring,2)
1578  if (natom/=natom_spring.or.(nspring/=1.and.nspring/=trotter)) then
1579    msg='Wrong dimensions for array spring !'
1580    MSG_BUG(msg)
1581  end if
1582 
1583  if (transform==0) then
1584    do iimage=1,trotter
1585      ispring=min(nspring,iimage)
1586      iimagep=iimage+1; iimagem=iimage-1
1587      if(iimage==trotter) iimagep=1
1588      if(iimage==1)       iimagem=trotter
1589      do iatom=1,natom
1590        do ii=1,3
1591          forces(ii,iatom,iimage)= &
1592 &               forces(ii,iatom,iimage)/dble(trotter) &
1593 &             - spring(iatom,ispring)*(two*xcart(ii,iatom,iimage)-xcart(ii,iatom,iimagem) &
1594 &                                                                -xcart(ii,iatom,iimagep))
1595        end do
1596      end do
1597    end do
1598 
1599  else
1600    do iimage=1,trotter
1601      ispring=min(nspring,iimage)
1602      do iatom=1,natom
1603        do ii=1,3
1604          forces(ii,iatom,iimage)= &
1605 &               forces(ii,iatom,iimage)/dble(trotter) &
1606 &             - spring(iatom,ispring)*xcart(ii,iatom,iimage)
1607        end do
1608      end do
1609    end do
1610 
1611  end if
1612 
1613 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.

PARENTS

      gstateimg

CHILDREN

SOURCE

145 subroutine pimd_init(dtset,pimd_param,is_master)
146 
147 
148 !This section has been created automatically by the script Abilint (TD).
149 !Do not modify the following lines by hand.
150 #undef ABI_FUNC
151 #define ABI_FUNC 'pimd_init'
152 !End of the abilint section
153 
154  implicit none
155 
156 !Arguments ------------------------------------
157 !scalars
158  logical,intent(in) :: is_master
159  type(dataset_type),target,intent(in) :: dtset
160  type(pimd_type),intent(inout) :: pimd_param
161 !Local variables-------------------------------
162 !scalars
163  integer :: ierr
164  character(len=200) :: msg
165 
166 !************************************************************************
167 
168  call pimd_nullify(pimd_param)
169 
170  if((dtset%imgmov==9).or.(dtset%imgmov==10).or.(dtset%imgmov==13))then
171    pimd_param%adpimd      = dtset%adpimd
172    pimd_param%constraint  = dtset%pimd_constraint
173    pimd_param%irandom     = dtset%irandom
174    pimd_param%nnos        = dtset%nnos
175    pimd_param%ntypat      = dtset%ntypat
176    pimd_param%optcell     = dtset%optcell
177    pimd_param%pitransform = dtset%pitransform
178    pimd_param%adpimd_gamma= dtset%adpimd_gamma
179    pimd_param%vis         = dtset%vis
180    pimd_param%bmass       = dtset%bmass
181    pimd_param%dtion       = dtset%dtion
182    pimd_param%friction    = dtset%friction
183    pimd_param%mdtemp      =>dtset%mdtemp
184    pimd_param%pimass      =>dtset%pimass
185    pimd_param%strtarget   =>dtset%strtarget
186    pimd_param%amu         =>dtset%amu_orig(:,1)
187    pimd_param%qmass       =>dtset%qmass
188    pimd_param%typat       =>dtset%typat
189    pimd_param%wtatcon     =>dtset%wtatcon
190    if(dtset%imgmov==10)then
191      pimd_param%use_qtb=1
192      if(is_master)then
193        call pimd_init_qtb(dtset,pimd_param%qtb_file_unit)
194      end if
195    end if
196    if(dtset%imgmov==13)then
197      ABI_ALLOCATE(pimd_param%zeta_prev,(3,dtset%natom,dtset%nimage,dtset%nnos))
198      ABI_ALLOCATE(pimd_param%zeta     ,(3,dtset%natom,dtset%nimage,dtset%nnos))
199      ABI_ALLOCATE(pimd_param%zeta_next,(3,dtset%natom,dtset%nimage,dtset%nnos))
200      ABI_ALLOCATE(pimd_param%dzeta    ,(3,dtset%natom,dtset%nimage,dtset%nnos))
201      pimd_param%zeta_prev=zero
202      pimd_param%zeta     =zero
203      pimd_param%zeta_next=zero
204      pimd_param%dzeta    =zero
205    end if
206    if(dtset%useria==37)then
207      ierr=open_file('pimd_traj.dat',msg,newunit=pimd_param%traj_unit,form='unformatted')
208      if (ierr/=0) then
209        MSG_ERROR(msg)
210      end if
211    end if
212  end if
213 
214 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.

PARENTS

      m_pimd

CHILDREN

SOURCE

386 subroutine pimd_init_qtb(dtset,qtb_file_unit)
387 
388 
389 !This section has been created automatically by the script Abilint (TD).
390 !Do not modify the following lines by hand.
391 #undef ABI_FUNC
392 #define ABI_FUNC 'pimd_init_qtb'
393 !End of the abilint section
394 
395  implicit none
396 
397 !Arguments ------------------------------------
398 !scalars
399  integer,intent(out) :: qtb_file_unit
400  type(dataset_type),target,intent(in) :: dtset
401 !Local variables-------------------------------
402 !scalars
403  integer :: ierr,ndof_qtb,ntimimage_qtb,nimage_qtb
404  real(dp) :: dtion_qtb,mdtemp_qtb
405  character(len=200) :: msg
406 
407 !************************************************************************
408 
409 !Try to open PIQTB random force file
410  ierr=open_file('piqtb_force',msg,newunit=qtb_file_unit,&
411 &               form='unformatted',status='old')
412 
413 !Read first line of the file
414  read(qtb_file_unit) dtion_qtb,ntimimage_qtb,mdtemp_qtb,nimage_qtb,ndof_qtb
415 
416 !Check consistency of the read parameters with ABINIT input file
417  if (abs(dtion_qtb-dtset%dtion)>tol6) then
418    msg='dtion read from piqtb_force file different from dtion in input file!'
419    MSG_ERROR(msg)
420  end if
421  if (abs(mdtemp_qtb-dtset%mdtemp(2))>tol6) then
422    msg='mdtemp read from piqtb_force file different from mdtemp(2) in input file!'
423    MSG_ERROR(msg)
424  end if
425  if (ntimimage_qtb<dtset%ntimimage) then
426    msg='ntimimage read from piqtb_force file smaller than ntimimage in input file!'
427    MSG_ERROR(msg)
428  end if
429  if (nimage_qtb/=dtset%nimage) then
430    msg='nimage read from piqtb_force file different from nimage in input file!'
431    MSG_ERROR(msg)
432  end if
433  if (ndof_qtb/=3*dtset%natom*dtset%nimage) then
434    msg='Nb of degrees of freedom read from piqtb_force not consistent with input file!'
435    MSG_ERROR(msg)
436  end if
437 
438 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

PARENTS

      pimd_langevin_npt,pimd_langevin_nvt,pimd_nosehoover_npt
      pimd_nosehoover_nvt

CHILDREN

SOURCE

 955 subroutine pimd_initvel(iseed,mass,natom,temperature,trotter,vel,constraint,wtatcon)
 956 
 957 
 958 !This section has been created automatically by the script Abilint (TD).
 959 !Do not modify the following lines by hand.
 960 #undef ABI_FUNC
 961 #define ABI_FUNC 'pimd_initvel'
 962 !End of the abilint section
 963 
 964  implicit none
 965 
 966 !Arguments ------------------------------------
 967 !scalars
 968  integer,intent(in) :: constraint,natom,trotter
 969  integer,intent(inout) :: iseed
 970  real(dp),intent(in) :: temperature
 971 !arrays
 972  real(dp),intent(in) :: mass(:,:),wtatcon(3,natom)
 973  real(dp),intent(out) :: vel(3,natom,trotter)
 974 !Local variables-------------------------------
 975 !scalars
 976  integer :: iatom,ii,iimage,imass,natom_mass,nmass
 977  real(dp) :: mtot,rescale_vel
 978  character(len=500) :: msg
 979 !arrays
 980  real(dp) :: mvini(3)
 981 
 982 !************************************************************************
 983 
 984  natom_mass=size(mass,1);nmass=size(mass,2)
 985  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
 986    msg='Wrong dimensions for array mass !'
 987    MSG_BUG(msg)
 988  end if
 989 
 990 !Compute total mass (of non-constrained atoms)
 991  if (constraint==0) then
 992    mtot=sum(mass(1:natom,1:nmass))
 993  else
 994    mtot=zero
 995    do iatom=1,natom
 996      if (all(abs(wtatcon(:,iatom))<tol8)) mtot=mtot+sum(mass(iatom,1:nmass))
 997    end do
 998  end if
 999  if (nmass==1) mtot=mtot*dble(trotter)
1000 
1001 !Initialize randomly the velocities
1002  do iimage=1,trotter
1003    imass=min(nmass,iimage)
1004    do iatom=1,natom
1005      do ii=1,3
1006        vel(ii,iatom,iimage)=sqrt(kb_HaK*temperature/mass(iatom,imass))*cos(two_pi*uniformrandom(iseed))
1007        vel(ii,iatom,iimage)=vel(ii,iatom,iimage)*sqrt(-two*log(uniformrandom(iseed)))
1008      end do
1009    end do
1010  end do
1011 
1012 !Cancel velocities of constrained atoms
1013  if (constraint/=0) then
1014    do iimage=1,trotter
1015      do iatom=1,natom
1016        if (any(abs(wtatcon(:,iatom))>=tol8)) vel(:,iatom,iimage)=zero
1017      end do
1018    end do
1019  end if
1020 
1021 !Make sure that the (sum of m_i v_i) at step zero is zero
1022  mvini=zero
1023  do iimage=1,trotter
1024    imass=min(nmass,iimage)
1025    do iatom=1,natom
1026      do ii=1,3
1027       mvini(ii)=mvini(ii)+mass(iatom,imass)*vel(ii,iatom,iimage)
1028      end do
1029    end do
1030  end do
1031  if (constraint==0) then
1032    do iimage=1,trotter
1033      do iatom=1,natom
1034        do ii=1,3
1035          vel(ii,iatom,iimage)=vel(ii,iatom,iimage)-(mvini(ii)/mtot)
1036        end do
1037      end do
1038    end do
1039  else
1040    do iimage=1,trotter
1041      do iatom=1,natom
1042        if (all(abs(wtatcon(:,iatom))<tol8)) vel(:,iatom,iimage)=vel(:,iatom,iimage)-(mvini(:)/mtot)
1043      end do
1044    end do
1045  end if
1046 
1047 !Now rescale the velocities to give the exact temperature
1048  rescale_vel=sqrt(temperature/pimd_temperature(mass,vel))
1049  vel(:,:,:)=vel(:,:,:)*rescale_vel
1050 
1051 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

PARENTS

CHILDREN

SOURCE

524 function pimd_is_restart(mass,vel,vel_cell)
525 
526 
527 !This section has been created automatically by the script Abilint (TD).
528 !Do not modify the following lines by hand.
529 #undef ABI_FUNC
530 #define ABI_FUNC 'pimd_is_restart'
531 !End of the abilint section
532 
533  implicit none
534 
535 !Arguments ------------------------------------
536 !scalars
537  integer :: pimd_is_restart
538 !arrays
539  real(dp),intent(in) :: mass(:,:),vel(:,:,:)
540  real(dp),intent(in),optional :: vel_cell(:,:,:)
541 !Local variables-------------------------------
542 !scalars
543  real(dp),parameter :: zero_temp=tol7
544 !arrays
545 
546 !************************************************************************
547 
548  pimd_is_restart=0
549  if (pimd_temperature(mass,vel)>zero_temp) pimd_is_restart=1
550  if (present(vel_cell)) then
551    if (maxval(vel_cell)>zero_temp) pimd_is_restart=pimd_is_restart+10
552  end if
553 
554 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

PARENTS

      pimd_langevin_npt,pimd_langevin_nvt

CHILDREN

SOURCE

1647 subroutine pimd_langevin_forces(alea,forces,forces_langevin,friction,&
1648 &                               langev,mass,natom,trotter,vel)
1649 
1650 
1651 !This section has been created automatically by the script Abilint (TD).
1652 !Do not modify the following lines by hand.
1653 #undef ABI_FUNC
1654 #define ABI_FUNC 'pimd_langevin_forces'
1655 !End of the abilint section
1656 
1657  implicit none
1658 
1659 !Arguments ------------------------------------
1660 !scalars
1661  integer,intent(in) :: natom,trotter
1662  real(dp),intent(in) :: friction
1663 !arrays
1664  real(dp),intent(in) :: alea(3,natom,trotter),forces(3,natom,trotter)
1665  real(dp),intent(in) :: langev(:,:),mass(:,:),vel(3,natom,trotter)
1666  real(dp),intent(out) :: forces_langevin(3,natom,trotter)
1667 !Local variables-------------------------------
1668 !scalars
1669  integer :: iatom,ii,iimage,imass,natom_mass,nmass
1670  character(len=500) :: msg
1671 !arrays
1672 
1673 !************************************************************************
1674 
1675  natom_mass=size(mass,1);nmass=size(mass,2)
1676  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
1677    msg='Wrong dimensions for array mass !'
1678    MSG_BUG(msg)
1679  end if
1680 
1681  do iimage=1,trotter
1682    imass=min(nmass,iimage)
1683    do iatom=1,natom
1684      do ii=1,3
1685        forces_langevin(ii,iatom,iimage)=forces(ii,iatom,iimage) &
1686 &                    + langev(iatom,imass)*alea(ii,iatom,iimage) &
1687 &                    - friction*mass(iatom,imass)*vel(ii,iatom,iimage)
1688      end do
1689    end do
1690  end do
1691 
1692 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)

PARENTS

      pimd_langevin_npt,pimd_langevin_nvt

CHILDREN

SOURCE

1087 subroutine pimd_langevin_random(alea,irandom,iseed,langev,mass,natom,trotter,zeroforce)
1088 
1089 
1090 !This section has been created automatically by the script Abilint (TD).
1091 !Do not modify the following lines by hand.
1092 #undef ABI_FUNC
1093 #define ABI_FUNC 'pimd_langevin_random'
1094 !End of the abilint section
1095 
1096  implicit none
1097 
1098 !Arguments ------------------------------------
1099 !scalars
1100  integer,intent(in) :: irandom,natom,trotter,zeroforce
1101  integer,intent(inout) :: iseed
1102 !arrays
1103  real(dp),intent(in) :: langev(:,:),mass(:,:)
1104  real(dp),intent(out) :: alea(3,natom,trotter)
1105 !Local variables-------------------------------
1106 !scalars
1107  integer :: iatom,ii,iimage,imass,nmass,natom_mass
1108  real(dp) :: mtot,r1,r2
1109  character(len=500) :: msg
1110 !arrays
1111  real(dp) :: total(3)
1112 
1113 !************************************************************************
1114 
1115  natom_mass=size(mass,1);nmass=size(mass,2)
1116  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
1117    msg='Wrong dimensions for array mass !'
1118    MSG_BUG(msg)
1119  end if
1120 
1121  mtot=sum(mass(1:natom,1:nmass))
1122  if (nmass==1) mtot=mtot*trotter
1123 
1124 !Draw random numbers
1125  do iimage=1,trotter
1126    do iatom=1,natom
1127      do ii=1,3
1128        select case(irandom)
1129        case(1)
1130          r1=uniformrandom(iseed)
1131          r2=uniformrandom(iseed)
1132        case(2)
1133          call random_number(r1)
1134          call random_number(r2)
1135        case(3)
1136          r1=ZBQLU01(zero)
1137          r2=ZBQLU01(zero)
1138        end select
1139        alea(ii,iatom,iimage)= cos(two_pi*r1)*sqrt(-log(r2)*two)
1140      end do
1141    end do
1142  end do
1143 
1144 !Make sure that the sum of random forces is zero
1145  if(zeroforce==1)then
1146    total=zero
1147    mtot=zero
1148    do iimage=1,trotter
1149      imass=min(nmass,iimage)
1150      do iatom=1,natom
1151        mtot=mtot+mass(iatom,imass)
1152      end do
1153    end do
1154    do iimage=1,trotter
1155      imass=min(nmass,iimage)
1156      do iatom=1,natom
1157        do ii=1,3
1158          total(ii)=total(ii)+langev(iatom,imass)*alea(ii,iatom,iimage)
1159        end do
1160      end do
1161    end do
1162    do iimage=1,trotter
1163      imass=min(nmass,iimage)
1164      do iatom=1,natom
1165        do ii=1,3
1166          alea(ii,iatom,iimage)= alea(ii,iatom,iimage)- &
1167 &        (total(ii)*mass(iatom,imass))/(langev(iatom,imass)*mtot)
1168        end do
1169      end do
1170    end do
1171 !  now random forces have been rescaled so that their sum is zero
1172  end if
1173 
1174 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)

PARENTS

      pimd_langevin_npt

CHILDREN

SOURCE

1307 subroutine pimd_langevin_random_bar(alea_bar,irandom,iseed)
1308 
1309 
1310 !This section has been created automatically by the script Abilint (TD).
1311 !Do not modify the following lines by hand.
1312 #undef ABI_FUNC
1313 #define ABI_FUNC 'pimd_langevin_random_bar'
1314 !End of the abilint section
1315 
1316  implicit none
1317 
1318 !Arguments ------------------------------------
1319 !scalars
1320  integer,intent(in) :: irandom
1321  integer,intent(inout) :: iseed
1322 !arrays
1323  real(dp),intent(out) :: alea_bar(3,3)
1324 !Local variables-------------------------------
1325 !scalars
1326  integer :: ii,jj
1327  real(dp) :: r1,r2
1328 !arrays
1329 
1330 !************************************************************************
1331 
1332 !Draw random numbers
1333  do ii=1,3
1334    do jj=1,3
1335      select case(irandom)
1336      case(1)
1337        r1=uniformrandom(iseed)
1338        r2=uniformrandom(iseed)
1339      case(2)
1340        call random_number(r1)
1341        call random_number(r2)
1342      case(3)
1343        r1=ZBQLU01(zero)
1344        r2=ZBQLU01(zero)
1345      end select
1346      alea_bar(ii,jj)= cos(two_pi*r1)*sqrt(-log(r2)*two)
1347    end do
1348  end do
1349 
1350 !Symmetrize
1351  alea_bar(1,2)=half*(alea_bar(1,2)+alea_bar(2,1))
1352  alea_bar(1,3)=half*(alea_bar(1,3)+alea_bar(3,1))
1353  alea_bar(2,3)=half*(alea_bar(2,3)+alea_bar(3,2))
1354  alea_bar(2,1)=alea_bar(1,2)
1355  alea_bar(3,1)=alea_bar(1,3)
1356  alea_bar(3,2)=alea_bar(2,3)
1357 
1358 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)

PARENTS

      pimd_langevin_npt,pimd_langevin_nvt

CHILDREN

SOURCE

1388 subroutine pimd_langevin_random_init(irandom,iseed)
1389 
1390 
1391 !This section has been created automatically by the script Abilint (TD).
1392 !Do not modify the following lines by hand.
1393 #undef ABI_FUNC
1394 #define ABI_FUNC 'pimd_langevin_random_init'
1395 !End of the abilint section
1396 
1397  implicit none
1398 
1399 !Arguments ------------------------------------
1400 !scalars
1401  integer,intent(in) :: irandom
1402  integer,intent(inout) :: iseed
1403 
1404 !************************************************************************
1405 
1406  if (irandom==3) then
1407    call ZBQLINI(0)
1408  end if
1409 
1410 !Fake statement
1411  return;if (.false.) iseed=zero
1412 
1413 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

PARENTS

      pimd_langevin_nvt

CHILDREN

SOURCE

1204 subroutine pimd_langevin_random_qtb(alea,langev,mass,natom,qtb_file_unit,trotter,zeroforce)
1205 
1206 
1207 !This section has been created automatically by the script Abilint (TD).
1208 !Do not modify the following lines by hand.
1209 #undef ABI_FUNC
1210 #define ABI_FUNC 'pimd_langevin_random_qtb'
1211 !End of the abilint section
1212 
1213  implicit none
1214 
1215 !Arguments ------------------------------------
1216 !scalars
1217  integer,intent(in) :: natom,qtb_file_unit,trotter,zeroforce
1218 !arrays
1219  real(dp),intent(in) :: langev(:,:),mass(:,:)
1220  real(dp),intent(out) :: alea(3,natom,trotter)
1221 !Local variables-------------------------------
1222 !scalars
1223  integer :: iatom,ii,iimage,imass,nmass,natom_mass
1224  real(dp) :: mtot
1225  character(len=500) :: msg
1226 !arrays
1227  real(sp) :: alea_sp(3,natom,trotter)
1228  real(dp) :: total(3)
1229 
1230 !************************************************************************
1231 
1232  if (qtb_file_unit<0) then
1233    msg='QTB forces file unit should be positive!'
1234    MSG_BUG(msg)
1235  end if
1236  natom_mass=size(mass,1);nmass=size(mass,2)
1237  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
1238    msg='Wrong dimensions for array mass !'
1239    MSG_BUG(msg)
1240  end if
1241 
1242 !Read QTB random forces
1243  read(qtb_file_unit) alea_sp(1:3,1:natom,1:trotter)
1244  alea(:,:,:)=dble(alea_sp(:,:,:))
1245 
1246 !Make sure that the sum of random forces is zero
1247  if(zeroforce==1)then
1248    total=zero
1249    mtot=zero
1250    do iimage=1,trotter
1251      imass=min(nmass,iimage)
1252      do iatom=1,natom
1253        mtot=mtot+mass(iatom,imass)
1254      end do
1255    end do
1256    do iimage=1,trotter
1257      imass=min(nmass,iimage)
1258      do iatom=1,natom
1259        do ii=1,3
1260          total(ii)=total(ii)+langev(iatom,imass)*alea(ii,iatom,iimage)
1261        end do
1262      end do
1263    end do
1264    do iimage=1,trotter
1265      imass=min(nmass,iimage)
1266      do iatom=1,natom
1267        do ii=1,3
1268          alea(ii,iatom,iimage)= alea(ii,iatom,iimage)- &
1269 &        (total(ii)*mass(iatom,imass))/(langev(iatom,imass)*mtot)
1270        end do
1271      end do
1272    end do
1273 !  now random forces have been rescaled so that their sum is zero
1274  end if
1275 
1276 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 !

PARENTS

      pimd_langevin_npt,pimd_langevin_nvt,pimd_nosehoover_nvt

CHILDREN

SOURCE

2849 subroutine pimd_mass_spring(inertmass,kt,mass,natom,quantummass,spring,transform,trotter)
2850 
2851 
2852 !This section has been created automatically by the script Abilint (TD).
2853 !Do not modify the following lines by hand.
2854 #undef ABI_FUNC
2855 #define ABI_FUNC 'pimd_mass_spring'
2856 !End of the abilint section
2857 
2858  implicit none
2859 
2860 !Arguments ------------------------------------
2861 !scalars
2862  integer,intent(in) :: natom,transform,trotter
2863  real(dp),intent(in) :: kt
2864 !arrays
2865  real(dp),intent(in) :: inertmass(natom),quantummass(natom)
2866  real(dp),intent(out) :: mass(:,:),spring(:,:)
2867 !Local variables-------------------------------
2868 !scalars
2869  integer :: iimage,kk,natom_mass,natom_spring,nmass,nspring
2870  real(dp) :: gammasquare
2871  character(len=500) :: msg
2872 !arrays
2873  real(dp),allocatable :: mass_temp(:,:),lambda(:)
2874 
2875 !************************************************************************
2876 
2877  natom_mass  =size(mass  ,1);nmass  =size(mass  ,2)
2878  natom_spring=size(spring,1);nspring=size(spring,2)
2879  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
2880    msg='Wrong dimensions for array mass !'
2881    MSG_BUG(msg)
2882  end if
2883  if (natom/=natom_spring.or.(nspring/=1.and.nspring/=trotter)) then
2884    msg='Wrong dimensions for array spring !'
2885    MSG_BUG(msg)
2886  end if
2887 
2888 !=== No transformation ===================================================
2889  if (transform==0) then
2890    !2nd dim of mass and spring = 1
2891    mass(1:natom,1)=inertmass(1:natom)
2892    spring(1:natom,1)=quantummass(1:natom)*dble(trotter)*kt*kt
2893 
2894 !=== Normal mode transformation ==========================================
2895  else if (transform==1) then
2896 
2897    ABI_ALLOCATE(lambda,(trotter))
2898    lambda(1)=zero; lambda(trotter)=four*dble(trotter)
2899    do kk=2,trotter/2
2900      lambda(2*kk-2)=four*dble(trotter)* &
2901 &                  (one-cos(two_pi*dble(kk-1)/dble(trotter)))
2902      lambda(2*kk-1)=four*dble(trotter)* &
2903 &                  (one-cos(two_pi*dble(kk-1)/dble(trotter)))
2904    end do
2905 
2906    !normal mode masses
2907    do iimage=1,trotter
2908      mass(:,iimage)=quantummass(:)*lambda(iimage)
2909    end do
2910 
2911    do iimage=1,trotter
2912      spring(:,iimage)=mass(:,iimage)*dble(trotter)*kt*kt
2913    end do
2914 
2915    !fictitious masses
2916    mass(:,1)=inertmass(:)
2917 
2918    !from 2 to P not changed except if adiabatic PIMD
2919    !see Hone et al, JCP 124, 154103 (2006)
2920    gammasquare=one  !adiabaticity parameter
2921    do iimage=2,trotter
2922      mass(:,iimage)=mass(:,iimage)/gammasquare
2923    end do
2924 
2925    ABI_DEALLOCATE(lambda)
2926 
2927 !=== Staging transformation ==============================================
2928  else if (transform==2) then
2929 
2930    !Fictitious masses
2931    mass(1:natom,1)=inertmass(1:natom)
2932    if (nmass>1) then
2933      do iimage=2,trotter
2934        mass(1:natom,iimage)=inertmass(1:natom)*dble(iimage)/dble(iimage-1)
2935      end do
2936    end if
2937 
2938    !Staging masses (mass_temp)
2939    ABI_ALLOCATE(mass_temp,(natom,trotter))
2940    mass_temp(1:natom,1)=zero
2941    if (nmass>1) then
2942      do iimage=2,trotter
2943        mass_temp(1:natom,iimage)=quantummass(1:natom)*dble(iimage)/dble(iimage-1)
2944      end do
2945    end if
2946 
2947    spring(1:natom,1)=mass_temp(1:natom,1)*dble(trotter)*kt*kt
2948    if (nspring>1) then
2949      do iimage=2,trotter
2950        spring(1:natom,iimage)=mass_temp(1:natom,iimage)*dble(trotter)*kt*kt
2951      end do
2952    end if
2953    ABI_DEALLOCATE(mass_temp)
2954 
2955  end if
2956 
2957 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

PARENTS

      pimd_nosehoover_nvt

CHILDREN

SOURCE

1727 subroutine pimd_nosehoover_forces(dzeta,forces,forces_nosehoover,mass,natom,&
1728 &                                 nnos,trotter,vel)
1729 
1730 
1731 !This section has been created automatically by the script Abilint (TD).
1732 !Do not modify the following lines by hand.
1733 #undef ABI_FUNC
1734 #define ABI_FUNC 'pimd_nosehoover_forces'
1735 !End of the abilint section
1736 
1737  implicit none
1738 
1739 !Arguments ------------------------------------
1740 !scalars
1741  integer,intent(in) :: natom,nnos,trotter
1742 !arrays
1743  real(dp),intent(in) :: dzeta(3,natom,trotter,nnos),forces(3,natom,trotter)
1744  real(dp),intent(in) :: vel(3,natom,trotter)
1745  real(dp),intent(in) :: mass(:,:)
1746  real(dp),intent(out) :: forces_nosehoover(3,natom,trotter)
1747 !Local variables-------------------------------
1748 !scalars
1749  integer :: iatom,ii,iimage,imass,natom_mass,nmass
1750  character(len=500) :: msg
1751 !arrays
1752 
1753 !************************************************************************
1754 
1755  natom_mass=size(mass,1);nmass=size(mass,2)
1756  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
1757    msg='Wrong dimensions for array mass !'
1758    MSG_BUG(msg)
1759  end if
1760 
1761   do iimage=1,trotter
1762     imass=min(nmass,iimage)
1763     do iatom=1,natom
1764       do ii=1,3
1765         forces_nosehoover(ii,iatom,iimage)=forces(ii,iatom,iimage) &
1766 &           - mass(iatom,imass)*dzeta(ii,iatom,iimage,1)*vel(ii,iatom,iimage)
1767       end do
1768     end do
1769   end do
1770 
1771 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

PARENTS

      pimd_nosehoover_nvt

CHILDREN

SOURCE

2196 subroutine pimd_nosehoover_propagate(dtion,dzeta,mass,natom,nnos,qmass,temperature,&
2197 &                                    trotter,vel,zeta,zeta_next,zeta_prev,itimimage,transform)
2198 
2199 
2200 !This section has been created automatically by the script Abilint (TD).
2201 !Do not modify the following lines by hand.
2202 #undef ABI_FUNC
2203 #define ABI_FUNC 'pimd_nosehoover_propagate'
2204 !End of the abilint section
2205 
2206  implicit none
2207 
2208 !Arguments ------------------------------------
2209 !scalars
2210  integer,intent(in) :: itimimage,natom,nnos,transform,trotter
2211  real(dp),intent(in) :: dtion,temperature
2212 !arrays
2213  real(dp),intent(in) :: qmass(nnos),vel(3,natom,trotter)
2214  real(dp),intent(in) :: mass(:,:)
2215  real(dp),intent(in) :: dzeta(3,natom,trotter,nnos),zeta_prev(3,natom,trotter,nnos)
2216  real(dp),intent(in) :: zeta(3,natom,trotter,nnos)
2217  real(dp),intent(out) :: zeta_next(3,natom,trotter,nnos)
2218 !Local variables-------------------------------
2219 !scalars
2220  integer :: iatom,ii,iimage,inos,natom_mass,nmass
2221  real(dp) :: kt
2222  character(len=500) :: msg
2223 !arrays
2224  real(dp),allocatable :: thermforces(:,:,:,:)
2225 
2226 !************************************************************************
2227 
2228  natom_mass=size(mass,1);nmass=size(mass,2)
2229  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
2230    msg='Wrong dimensions for array mass !'
2231    MSG_BUG(msg)
2232  end if
2233  if (nnos<3) then
2234    msg='Not available for nnos<3 !'
2235    MSG_BUG(msg)
2236  end if
2237 
2238  kt=temperature*kb_HaK
2239  ABI_ALLOCATE(thermforces,(3,natom,trotter,nnos))
2240 
2241 !forces on the thermostats
2242  do inos=1,nnos
2243    if(inos==1)then
2244      do iimage=1,trotter
2245        do iatom=1,natom
2246          do ii=1,3
2247            select case(transform)
2248            case(0)
2249              thermforces(ii,iatom,iimage,1)=mass(iatom,1)*(vel(ii,iatom,iimage)**2)-kt
2250            case(1)
2251              thermforces(ii,iatom,iimage,1)=mass(iatom,iimage)*(vel(ii,iatom,iimage)**2)-kt
2252            case(2)
2253              thermforces(ii,iatom,iimage,1)=mass(iatom,iimage)*(vel(ii,iatom,iimage)**2)-kt
2254            end select
2255          end do
2256        end do
2257      end do
2258    elseif(inos==(nnos-1))then
2259      do iimage=1,trotter
2260        do iatom=1,natom
2261          do ii=1,3
2262            thermforces(ii,iatom,iimage,inos)=&
2263 &          (qmass(nnos-2)*(dzeta(ii,iatom,iimage,nnos-2)**2)-kt) + &
2264 &          (qmass(nnos  )*(dzeta(ii,iatom,iimage,nnos)  **2)-kt)
2265          end do
2266        end do
2267      end do
2268    elseif(inos==nnos)then
2269      do iimage=1,trotter
2270        do iatom=1,natom
2271          do ii=1,3
2272            thermforces(ii,iatom,iimage,inos)=&
2273 &          qmass(nnos-1)*(dzeta(ii,iatom,iimage,nnos-1)**2)-kt
2274          end do
2275        end do
2276      end do
2277    else
2278      do iimage=1,trotter
2279        do iatom=1,natom
2280          do ii=1,3
2281            thermforces(ii,iatom,iimage,inos)=&
2282 &          qmass(inos-1)*(dzeta(ii,iatom,iimage,inos-1)**2)-kt
2283          end do
2284        end do
2285      end do
2286    end if
2287  end do
2288 
2289  select case(itimimage)
2290  case(1) !taylor
2291 
2292  do inos=1,nnos
2293    if(inos==1)then
2294      zeta_next(:,:,:,1)=zeta(:,:,:,1)+ dzeta(:,:,:,1)*dtion + &
2295 &       (thermforces(:,:,:,1)-qmass(1)*dzeta(:,:,:,1)*dzeta(:,:,:,2))* &
2296 &       dtion*dtion/(two*qmass(1))
2297    elseif(inos==(nnos-1))then
2298      zeta_next(:,:,:,inos)=zeta(:,:,:,inos)+ dzeta(:,:,:,inos)*dtion + &
2299 &       (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,inos)*dzeta(:,:,:,nnos))* &
2300 &       dtion*dtion/(two*qmass(inos))
2301    elseif(inos==nnos)then
2302      zeta_next(:,:,:,inos)=zeta(:,:,:,inos)+ dzeta(:,:,:,inos)*dtion + &
2303 &       (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,nnos-1)*dzeta(:,:,:,nnos))* &
2304 &       dtion*dtion/(two*qmass(inos))
2305    else
2306      zeta_next(:,:,:,inos)=zeta(:,:,:,inos)+ dzeta(:,:,:,inos)*dtion + &
2307 &       (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,inos)*dzeta(:,:,:,inos+1))* &
2308 &       dtion*dtion/(two*qmass(inos))
2309    end if
2310  end do
2311 
2312  case default !verlet
2313 
2314  do inos=1,nnos
2315    if(inos==1)then
2316      zeta_next(:,:,:,1)=two*zeta(:,:,:,1) - zeta_prev(:,:,:,1) + &
2317 &       (thermforces(:,:,:,1)-qmass(1)*dzeta(:,:,:,1)*dzeta(:,:,:,2))* &
2318 &       dtion*dtion/qmass(1)
2319    elseif(inos==(nnos-1))then
2320      zeta_next(:,:,:,inos)=two*zeta(:,:,:,inos) - zeta_prev(:,:,:,inos) + &
2321 &       (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,inos)*dzeta(:,:,:,nnos))* &
2322 &       dtion*dtion/qmass(inos)
2323    elseif(inos==nnos)then
2324      zeta_next(:,:,:,inos)=two*zeta(:,:,:,inos) - zeta_prev(:,:,:,inos) + &
2325 &       (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,nnos-1)*dzeta(:,:,:,nnos))* &
2326 &       dtion*dtion/qmass(inos)
2327    else
2328      zeta_next(:,:,:,inos)=two*zeta(:,:,:,inos) - zeta_prev(:,:,:,inos) + &
2329 &       (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,inos)*dzeta(:,:,:,inos+1))* &
2330 &       dtion*dtion/qmass(inos)
2331    end if
2332  end do
2333 
2334  end select
2335 
2336  ABI_DEALLOCATE(thermforces)
2337 
2338 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.

PARENTS

      m_pimd

CHILDREN

SOURCE

241 subroutine pimd_nullify(pimd_param)
242 
243 
244 !This section has been created automatically by the script Abilint (TD).
245 !Do not modify the following lines by hand.
246 #undef ABI_FUNC
247 #define ABI_FUNC 'pimd_nullify'
248 !End of the abilint section
249 
250  implicit none
251 
252 !Arguments ------------------------------------
253 !scalars
254  type(pimd_type),intent(inout) :: pimd_param
255 
256 !************************************************************************
257 
258  pimd_param%adpimd       =  0
259  pimd_param%constraint   =  0
260  pimd_param%irandom      = -1
261  pimd_param%nnos         = -1
262  pimd_param%ntypat       = -1
263  pimd_param%optcell      = -1
264  pimd_param%pitransform  = -1
265  pimd_param%qtb_file_unit= -1
266  pimd_param%traj_unit    = -1
267  pimd_param%use_qtb      =  0
268  pimd_param%adpimd_gamma = one
269  pimd_param%vis          = zero
270  pimd_param%bmass        = zero
271  pimd_param%dtion        = zero
272  pimd_param%friction     = zero
273  nullify(pimd_param%mdtemp)
274  nullify(pimd_param%pimass)
275  nullify(pimd_param%strtarget)
276  nullify(pimd_param%amu)
277  nullify(pimd_param%qmass)
278  nullify(pimd_param%typat)
279  nullify(pimd_param%wtatcon)
280 
281 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

PARENTS

      pimd_langevin_nvt,pimd_nosehoover_nvt

CHILDREN

SOURCE

2029 subroutine pimd_predict_taylor(dtion,forces,mass,natom,trotter,vel,xcart,xcart_next)
2030 
2031 
2032 !This section has been created automatically by the script Abilint (TD).
2033 !Do not modify the following lines by hand.
2034 #undef ABI_FUNC
2035 #define ABI_FUNC 'pimd_predict_taylor'
2036 !End of the abilint section
2037 
2038  implicit none
2039 
2040 !Arguments ------------------------------------
2041 !scalars
2042  integer,intent(in) :: natom,trotter
2043  real(dp),intent(in) :: dtion
2044 !arrays
2045  real(dp),intent(in) :: forces(3,natom,trotter),vel(3,natom,trotter),xcart(3,natom,trotter)
2046  real(dp),intent(in) :: mass(:,:)
2047  real(dp),intent(out) :: xcart_next(3,natom,trotter)
2048 !Local variables-------------------------------
2049 !scalars
2050  integer :: iatom,ii,iimage,imass,natom_mass,nmass
2051  character(len=500) :: msg
2052 !arrays
2053 
2054 !************************************************************************
2055 
2056  natom_mass=size(mass,1);nmass=size(mass,2)
2057  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
2058    msg='Wrong dimensions for array mass !'
2059    MSG_BUG(msg)
2060  end if
2061 
2062  do iimage=1,trotter
2063    imass=min(nmass,iimage)
2064    do iatom=1,natom
2065      do ii=1,3
2066        xcart_next(ii,iatom,iimage)=xcart(ii,iatom,iimage) &
2067 &             + half*dtion*dtion*forces(ii,iatom,iimage)/mass(iatom,imass) &
2068 &             + dtion*vel(ii,iatom,iimage)
2069      end do
2070    end do
2071  end do
2072 
2073 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

PARENTS

      pimd_langevin_nvt,pimd_nosehoover_nvt

CHILDREN

SOURCE

2106 subroutine pimd_predict_verlet(dtion,forces,mass,natom,trotter,xcart,xcart_next,xcart_prev)
2107 
2108 
2109 !This section has been created automatically by the script Abilint (TD).
2110 !Do not modify the following lines by hand.
2111 #undef ABI_FUNC
2112 #define ABI_FUNC 'pimd_predict_verlet'
2113 !End of the abilint section
2114 
2115  implicit none
2116 
2117 !Arguments ------------------------------------
2118 !scalars
2119  integer,intent(in) :: natom,trotter
2120  real(dp),intent(in) :: dtion
2121 !arrays
2122  real(dp),intent(in) :: forces(3,natom,trotter)
2123  real(dp),intent(in) :: xcart(3,natom,trotter),xcart_prev(3,natom,trotter)
2124  real(dp),intent(in) :: mass(:,:)
2125  real(dp),intent(out) :: xcart_next(3,natom,trotter)
2126 !Local variables-------------------------------
2127 !scalars
2128  integer :: iatom,ii,iimage,imass,natom_mass,nmass
2129  character(len=500) :: msg
2130 !arrays
2131 
2132 !************************************************************************
2133 
2134  natom_mass=size(mass,1);nmass=size(mass,2)
2135  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
2136    msg='Wrong dimensions for array mass !'
2137    MSG_BUG(msg)
2138  end if
2139 
2140  do iimage=1,trotter
2141    imass=min(nmass,iimage)
2142    do iatom=1,natom
2143      do ii=1,3
2144        xcart_next(ii,iatom,iimage)= &
2145 &         two*xcart(ii,iatom,iimage) &
2146 &       - xcart_prev(ii,iatom,iimage) &
2147 &       + dtion*dtion*forces(ii,iatom,iimage)/mass(iatom,imass)
2148      end do
2149    end do
2150  end do
2151 
2152 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

PARENTS

      pimd_langevin_npt,pimd_langevin_nvt,pimd_nosehoover_npt
      pimd_nosehoover_nvt

CHILDREN

SOURCE

679 subroutine pimd_print(constraint,constraint_output,eharm,eharm_virial,epot,&
680 &          forces,inertmass,irestart,itimimage,kt,natom,optcell,prtstress,&
681 &          prtvolimg,rprimd,stress,temperature1,temperature2,traj_unit,&
682 &          trotter,vel,vel_cell,xcart,xred)
683 
684 
685 !This section has been created automatically by the script Abilint (TD).
686 !Do not modify the following lines by hand.
687 #undef ABI_FUNC
688 #define ABI_FUNC 'pimd_print'
689 !End of the abilint section
690 
691  implicit none
692 
693 !Arguments ------------------------------------
694 !scalars
695  integer,intent(in) :: constraint,irestart,itimimage,natom,optcell
696  integer,intent(in) :: prtstress,prtvolimg,traj_unit,trotter
697  real(dp),intent(in) :: eharm,eharm_virial,epot,kt,temperature1,temperature2
698 !arrays
699  real(dp),intent(in) :: constraint_output(2)
700  real(dp),intent(in) :: forces(3,natom,trotter),inertmass(natom)
701  real(dp),intent(in) :: rprimd(3,3),stress(3,3,3),vel(3,natom,trotter),vel_cell(3,3)
702  real(dp),intent(in) :: xcart(3,natom,trotter),xred(3,natom,trotter)
703 !Local variables-------------------------------
704 !scalars
705  integer :: iatom,ii,iimage
706  real(dp) :: mtot
707  character(len=500) :: msg
708 !arrays
709  real(dp) :: acell(3),cdm(3),forcetot(3),rprim(3,3)
710  real(dp),allocatable :: centroid(:,:),qudeloc(:)
711 
712 !************************************************************************
713 
714 !Temperature
715  if(itimimage==1)then
716    if(irestart==0) then
717      write(msg,'(2a)') ch10,' This is a PIMD calculation from scratch'
718    else if (irestart==1) then
719      write(msg,'(2a)') ch10,' This is a RESTART calculation'
720    end if
721    call wrtout(ab_out,msg,'COLL')
722    call wrtout(std_out,msg,'COLL')
723    write(msg,'(a,f12.5,a)') &
724 &   ' In the initial configuration, the temperature is ',temperature1,' K'
725    call wrtout(ab_out,msg,'COLL')
726    call wrtout(std_out,msg,'COLL')
727  end if
728  write(msg,'(2a,i5,a,f12.5,a)') ch10,&
729 &  ' At PIMD time step ',itimimage,', the temperature is',temperature2,' K'
730  call wrtout(ab_out,msg,'COLL')
731  call wrtout(std_out,msg,'COLL')
732 
733 !Energies
734  write(msg,'(4a,f18.9,a,a,a,f18.9,a,a)') ch10,&
735 &  ' Energy:',ch10, &
736 &  '   Internal energy (PRIMITIVE estimator) =',onehalf*dble(natom*trotter)*kt-eharm+epot ,' Ha',ch10, &
737 &  '   Internal energy (VIRIAL    estimator) =',onehalf*dble(natom)*kt+eharm_virial+epot,' Ha',ch10
738  call wrtout(ab_out,msg,'COLL')
739  call wrtout(std_out,msg,'COLL')
740 
741 !Stress tensor and pressure
742  write(msg,'(2a,3(2a,3f18.9))') ch10,&
743 &   ' Stress tensor from PRIMITIVE estimator (Ha/Bohr^3):',ch10, &
744 &   '   ',stress(1,1,1),stress(1,1,2),stress(1,1,3),ch10, &
745 &   '   ',stress(1,2,1),stress(1,2,2),stress(1,2,3),ch10, &
746 &   '   ',stress(1,3,1),stress(1,3,2),stress(1,3,3)
747  if (prtstress==1) then
748    call wrtout(ab_out,msg,'COLL')
749  end if
750  call wrtout(std_out,msg,'COLL')
751  write(msg,'(a,f18.9,a)') ' Pressure (primitive estimator) =', &
752 &  -third*(stress(1,1,1)+stress(1,2,2)+stress(1,3,3))*HaBohr3_GPa,' GPa'
753  call wrtout(ab_out,msg,'COLL')
754  call wrtout(std_out,msg,'COLL')
755 
756 !Data related to constraint eventually applied
757  if (constraint/=0) then
758    if (constraint==1) write(msg,'(2a)') ch10,' Blue Moon Ensemble method is activated:'
759    call wrtout(ab_out,msg,'COLL')
760    call wrtout(std_out,msg,'COLL')
761    write(msg,'(a,f18.10,2a,f18.10)') &
762 &     '  - Reaction coordinate =',constraint_output(1),ch10,&
763 &     '  - Instantaneous force on the reaction coord. =',constraint_output(2)
764    call wrtout(ab_out,msg,'COLL')
765    call wrtout(std_out,msg,'COLL')
766  end if
767 
768 !Total force
769  if (prtvolimg<=1) then
770    forcetot=zero
771    do iimage=1,trotter
772      do iatom=1,natom
773        do ii=1,3
774          forcetot(ii)=forcetot(ii)+forces(ii,iatom,iimage)
775        end do
776      end do
777    end do
778    write(msg,'(2a,3f18.10)') ch10,' Total force=',forcetot(1:3)
779    call wrtout(std_out,msg,'COLL')
780  end if
781 
782 !position of mass center
783  if (prtvolimg<=1) then
784    mtot=zero;cdm=zero
785    do iimage=1,trotter
786      do iatom=1,natom
787        cdm(:)=cdm(:)+inertmass(iatom)*xcart(:,iatom,iimage)
788        mtot=mtot+inertmass(iatom)
789      end do
790    end do
791    cdm=cdm/mtot
792    write(msg,'(3a,3f18.10)') ch10,' Center of mass:',ch10,cdm(:)
793    call wrtout(std_out,msg,'COLL')
794    call wrtout(ab_out,msg,'COLL')
795  end if
796 
797 !Positions
798  write(msg,'(2a)') ch10,' Atomic positions:'
799  call wrtout(std_out,msg,'COLL')
800  call wrtout(ab_out,msg,'COLL')
801  do iimage=1,trotter
802    select case(iimage)
803    case(1)
804     write(msg,'(a)') ' xred'
805    case(2,3,4,5,6,7,8,9)
806      write(msg,'(a,i1,a)') ' xred_',iimage,'img'
807    case(10:99)
808      write(msg,'(a,i2,a)') ' xred_',iimage,'img'
809    case default
810      write(msg,'(a,i3,a)') ' xred_',iimage,'img'
811    end select
812    call wrtout(std_out,msg,'COLL')
813    call wrtout(ab_out,msg,'COLL')
814    if (traj_unit>0) then
815      call wrtout(traj_unit,msg,'COLL')
816    end if
817    do iatom=1,natom
818      write(msg,'(3f18.10)') xred(1:3,iatom,iimage)
819      call wrtout(std_out,msg,'COLL')
820      call wrtout(ab_out,msg,'COLL')
821      if (traj_unit>0) then
822        call wrtout(traj_unit,msg,'COLL')
823      end if
824    end do
825  end do
826 
827 !Velocities
828  write(msg,'(2a)') ch10,' Velocities:'
829  call wrtout(std_out,msg,'COLL')
830  call wrtout(ab_out,msg,'COLL')
831  do iimage=1,trotter
832    select case(iimage)
833    case(1)
834      write(msg,'(a)') ' vel'
835    case(2,3,4,5,6,7,8,9)
836      write(msg,'(a,i1,a)') ' vel_',iimage,'img'
837    case(10:99)
838      write(msg,'(a,i2,a)') ' vel_',iimage,'img'
839    case default
840      write(msg,'(a,i3,a)') ' vel_',iimage,'img'
841    end select
842    call wrtout(std_out,msg,'COLL')
843    call wrtout(ab_out,msg,'COLL')
844    if (traj_unit>0) then
845      call wrtout(traj_unit,msg,'COLL')
846    end if
847    do iatom=1,natom
848      write(msg,'(3f18.10)') vel(1:3,iatom,iimage)
849      call wrtout(std_out,msg,'COLL')
850      call wrtout(ab_out,msg,'COLL')
851      if (traj_unit>0) then
852        call wrtout(traj_unit,msg,'COLL')
853      end if
854    end do
855  end do
856 
857  if (optcell>0) then
858 
859    call mkradim(acell,rprim,rprimd)
860 
861 !  Time derivative of rprimd
862    write(msg,'(2a)') ch10,' Time derivative of rprimd:'
863    call wrtout(std_out,msg,'COLL')
864    call wrtout(ab_out,msg,'COLL')
865    write(msg,'(2a,3(3f18.10))') ' vel_cell',ch10,vel_cell(:,:)
866    call wrtout(std_out,msg,'COLL')
867    call wrtout(ab_out,msg,'COLL')
868    if (traj_unit>0) then
869      call wrtout(traj_unit,msg,'COLL')
870    end if
871 
872 !  rprimd
873    write(msg,'(2a)') ch10,' Cell parameters:'
874    call wrtout(std_out,msg,'COLL')
875    call wrtout(ab_out,msg,'COLL')
876    write(msg,'(2a,3(3f18.10),3a,3f18.10)') ' rprim',ch10,rprim(:,:),ch10,&
877 &                                          ' acell',ch10,acell(:)
878    call wrtout(std_out,msg,'COLL')
879    call wrtout(ab_out,msg,'COLL')
880    if (traj_unit>0) then
881      call wrtout(traj_unit,msg,'COLL')
882    end if
883 
884  end if
885 
886 !Centroids and wave-packet spatial spreads
887  ABI_ALLOCATE(centroid,(3,natom))
888  ABI_ALLOCATE(qudeloc,(natom))
889  centroid=zero;qudeloc=zero
890  do iimage=1,trotter
891    do iatom=1,natom
892      do ii=1,3
893        centroid(ii,iatom)=centroid(ii,iatom)+xcart(ii,iatom,iimage)
894      end do
895    end do
896  end do
897  centroid=centroid/dble(trotter)
898  do iimage=1,trotter
899    do iatom=1,natom
900      do ii=1,3
901        qudeloc(iatom)=qudeloc(iatom)+((xcart(ii,iatom,iimage)-centroid(ii,iatom))**2)
902      end do
903    end do
904  end do
905  qudeloc(:)=sqrt(qudeloc(:)/dble(trotter))
906  write(msg,'(4a)') ch10,' Centroids and wave-packet spatial spreads (cart. coord.):',ch10,&
907 &  ' iat        centroid_x        centroid_y        centroid_z    spatial_spread'
908  call wrtout(std_out,msg,'COLL')
909  do iatom=1,natom
910    write(msg,'(i4,4f18.10)') iatom,centroid(1:3,iatom),qudeloc(iatom)
911    call wrtout(std_out,msg,'COLL')
912  end do
913  ABI_DEALLOCATE(centroid)
914  ABI_DEALLOCATE(qudeloc)
915 
916 !Fake statement
917  return;ii=prtvolimg
918 
919 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

PARENTS

      gstateimg

CHILDREN

SOURCE

464 subroutine pimd_skip_qtb(pimd_param)
465 
466 
467 !This section has been created automatically by the script Abilint (TD).
468 !Do not modify the following lines by hand.
469 #undef ABI_FUNC
470 #define ABI_FUNC 'pimd_skip_qtb'
471 !End of the abilint section
472 
473  implicit none
474 
475 !Arguments ------------------------------------
476 !scalars
477  type(pimd_type),intent(in) :: pimd_param
478 !arrays
479 !Local variables-------------------------------
480 !scalars
481  character(len=500) :: msg
482 !arrays
483 
484 !************************************************************************
485 
486  if (pimd_param%use_qtb==0) return
487 
488  if (pimd_param%qtb_file_unit<0) then
489    msg='QTB forces file unit should be positive!'
490    MSG_BUG(msg)
491  end if
492 
493 !Skip one line QTB random forces file
494  read(pimd_param%qtb_file_unit)
495 
496 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

PARENTS

      pimd_langevin_npt,pimd_langevin_nvt,pimd_nosehoover_npt
      pimd_nosehoover_nvt

CHILDREN

SOURCE

1810 subroutine pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,&
1811 &                        temperature,temperature_therm,trotter,vel,volume,xcart)
1812 
1813 
1814 !This section has been created automatically by the script Abilint (TD).
1815 !Do not modify the following lines by hand.
1816 #undef ABI_FUNC
1817 #define ABI_FUNC 'pimd_stresses'
1818 !End of the abilint section
1819 
1820  implicit none
1821 
1822 !Arguments ------------------------------------
1823 !scalars
1824  integer,intent(in) :: natom,trotter
1825  real(dp),intent(in) :: temperature,temperature_therm,volume
1826 !arrays
1827  real(dp),intent(in) :: vel(3,natom,trotter),xcart(3,natom,trotter)
1828  real(dp),intent(in) :: mass(:,:),stressin(3,3,trotter),quantummass(natom)
1829  real(dp),intent(out) :: stress_pimd(3,3,3)
1830 !Local variables-------------------------------
1831 !scalars
1832  integer :: iatom,ii,iimage,imass,jj,natom_mass,nmass,iimagep
1833  real(dp) :: stress_tmp(3,3,3),omega2,kt,kt_therm
1834  character(len=500) :: msg
1835 
1836 !************************************************************************
1837 
1838  natom_mass=size(mass,1);nmass=size(mass,2)
1839  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then
1840    msg='Wrong dimensions for array mass !'
1841    MSG_BUG(msg)
1842  end if
1843 
1844  kt=temperature*kb_HaK
1845  kt_therm=temperature_therm*kb_HaK
1846  stress_pimd=zero
1847 
1848 !I-PRIMITIVE ESTIMATOR
1849 !1-Kinetic part
1850  do ii=1,3
1851    stress_pimd(1,ii,ii)=dble(natom)*dble(trotter)*kt/volume
1852  end do
1853 
1854 !2-Potential part
1855  do iimage=1,trotter
1856    do ii=1,3
1857      do jj=1,3
1858        stress_pimd(1,ii,jj)=stress_pimd(1,ii,jj)-stressin(ii,jj,iimage)/dble(trotter)
1859        !minus to convert stress into pressure
1860      end do
1861    end do
1862  end do
1863 
1864 !3-Contribution from springs
1865  omega2=dble(trotter)*kt_therm*kt_therm
1866  do iimage=1,trotter
1867    iimagep=iimage+1
1868    if(iimage==trotter) then
1869      iimagep=1
1870    end if
1871    do iatom=1,natom
1872      do ii=1,3
1873        do jj=1,3
1874          stress_pimd(1,ii,jj)=stress_pimd(1,ii,jj)-quantummass(iatom)*omega2* &
1875 &        (xcart(ii,iatom,iimagep)-xcart(ii,iatom,iimage))*  &
1876 &        (xcart(jj,iatom,iimagep)-xcart(jj,iatom,iimage))/volume
1877        end do
1878      end do
1879    end do
1880  end do
1881 
1882 !II-Average of classical pressures
1883 !1-Kinetic part
1884  do iimage=1,trotter
1885   imass=min(nmass,iimage)
1886   do iatom=1,natom
1887     do ii=1,3
1888       do jj=1,3
1889         stress_pimd(2,ii,jj)=stress_pimd(2,ii,jj)+ &
1890 &          mass(iatom,imass)*vel(ii,iatom,iimage)*vel(jj,iatom,iimage)/volume
1891       end do
1892     end do
1893   end do
1894  end do
1895 
1896 !2-Contribution from electronic stress
1897  do iimage=1,trotter
1898    do ii=1,3
1899      do jj=1,3
1900        stress_pimd(2,ii,jj)=stress_pimd(2,ii,jj)-stressin(ii,jj,iimage)
1901      end do
1902    end do
1903  end do
1904  do ii=1,3
1905    do jj=1,3
1906      stress_pimd(2,ii,jj)=stress_pimd(2,ii,jj)/dble(trotter)
1907    end do
1908  end do
1909 
1910 !III-pressure from VIRIAL estimator: stress_pimd(3,:,:)
1911 !1-kinetic part
1912 ! do ii=1,3
1913 !   stress_pimd(3,ii,ii)=dfloat(natom)*kt/volume
1914 ! end do
1915 
1916 !Symmetrize internal pressure
1917  stress_tmp=stress_pimd
1918  stress_pimd(:,2,1)=half*(stress_tmp(:,1,2)+stress_tmp(:,2,1))
1919  stress_pimd(:,1,2)=half*(stress_tmp(:,1,2)+stress_tmp(:,2,1))
1920  stress_pimd(:,3,1)=half*(stress_tmp(:,1,3)+stress_tmp(:,3,1))
1921  stress_pimd(:,1,3)=half*(stress_tmp(:,1,3)+stress_tmp(:,3,1))
1922  stress_pimd(:,2,3)=half*(stress_tmp(:,3,2)+stress_tmp(:,2,3))
1923  stress_pimd(:,3,2)=half*(stress_tmp(:,3,2)+stress_tmp(:,2,3))
1924 
1925 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

PARENTS

CHILDREN

SOURCE

581 function pimd_temperature(mass,vel)
582 
583 
584 !This section has been created automatically by the script Abilint (TD).
585 !Do not modify the following lines by hand.
586 #undef ABI_FUNC
587 #define ABI_FUNC 'pimd_temperature'
588 !End of the abilint section
589 
590  implicit none
591 
592 !Arguments ------------------------------------
593 !scalars
594  real(dp) :: pimd_temperature
595 !arrays
596  real(dp),intent(in) :: mass(:,:),vel(:,:,:)
597 !Local variables-------------------------------
598 !scalars
599  integer :: iatom,idir,iimage,imass,natom,natom_mass,ndir,nimage,nmass
600  real(dp) :: v2
601  character(len=500) :: msg
602 !arrays
603 
604 !************************************************************************
605 
606  ndir=size(vel,1);natom=size(vel,2);nimage=size(vel,3)
607  natom_mass=size(mass,1);nmass=size(mass,2)
608  if (ndir/=3.or.natom<=0.or.nimage<=0) then
609    msg='Wrong sizes for vel array !'
610    MSG_BUG(msg)
611  end if
612  if (natom/=natom_mass.or.(nmass/=1.and.nmass/=nimage)) then
613    msg='Wrong dimensions for array mass !'
614    MSG_BUG(msg)
615  end if
616 
617  v2=zero
618  do iimage=1,nimage
619    imass=min(nmass,iimage)
620    do iatom=1,natom
621      do idir=1,3
622        v2=v2+vel(idir,iatom,iimage)*vel(idir,iatom,iimage)*mass(iatom,imass)
623      end do
624    end do
625  end do
626  pimd_temperature=v2/(dble(3*natom*nimage)*kb_HaK)
627 
628 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

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