Line data    Source code

       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2023 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief  Data type and methods dealing with PI calcs in normal mode coords
      10             : !> \author fawzi
      11             : !> \par    History
      12             : !>         2006-02 created
      13             : !>         2006-11 modified so it might actually work [hforbert]
      14             : !>         2009-04-07 moved from pint_types module to a separate file [lwalewski]
      15             : !>         2015-10 added alternative normal mode transformation needed by RPMD
      16             : !>                 [Felix Uhl
      17             : ! **************************************************************************************************
      18             : MODULE pint_normalmode
      19             :    USE input_constants,                 ONLY: propagator_cmd,&
      20             :                                               propagator_pimd,&
      21             :                                               propagator_rpmd
      22             :    USE input_section_types,             ONLY: section_vals_type,&
      23             :                                               section_vals_val_get
      24             :    USE kinds,                           ONLY: dp
      25             :    USE mathconstants,                   ONLY: pi,&
      26             :                                               twopi
      27             :    USE pint_types,                      ONLY: normalmode_env_type
      28             : #include "../base/base_uses.f90"
      29             : 
      30             :    IMPLICIT NONE
      31             :    PRIVATE
      32             : 
      33             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      34             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pint_normalmode'
      35             : 
      36             :    PUBLIC :: normalmode_env_create
      37             :    PUBLIC :: normalmode_release
      38             :    PUBLIC :: normalmode_init_masses
      39             :    PUBLIC :: normalmode_x2u
      40             :    PUBLIC :: normalmode_u2x
      41             :    PUBLIC :: normalmode_f2uf
      42             :    PUBLIC :: normalmode_calc_uf_h
      43             : 
      44             : CONTAINS
      45             : 
      46             : ! ***************************************************************************
      47             : !> \brief creates the data needed for a normal mode transformation
      48             : !> \param normalmode_env ...
      49             : !> \param normalmode_section ...
      50             : !> \param p ...
      51             : !> \param kT ...
      52             : !> \param propagator ...
      53             : !> \author Harald Forbert
      54             : ! **************************************************************************************************
      55          54 :    SUBROUTINE normalmode_env_create(normalmode_env, normalmode_section, p, kT, propagator)
      56             :       TYPE(normalmode_env_type), INTENT(OUT)             :: normalmode_env
      57             :       TYPE(section_vals_type), POINTER                   :: normalmode_section
      58             :       INTEGER, INTENT(in)                                :: p
      59             :       REAL(kind=dp), INTENT(in)                          :: kT
      60             :       INTEGER, INTENT(in)                                :: propagator
      61             : 
      62             :       INTEGER                                            :: i, j, k, li
      63             :       LOGICAL                                            :: explicit_gamma, explicit_modefactor
      64             :       REAL(kind=dp)                                      :: gamma_parameter, invsqrtp, pip, sqrt2p, &
      65             :                                                             twopip
      66             : 
      67         216 :       ALLOCATE (normalmode_env%x2u(p, p))
      68         216 :       ALLOCATE (normalmode_env%u2x(p, p))
      69         162 :       ALLOCATE (normalmode_env%lambda(p))
      70             : 
      71          54 :       normalmode_env%p = p
      72             : 
      73             :       CALL section_vals_val_get(normalmode_section, "Q_CENTROID", &
      74          54 :                                 r_val=normalmode_env%Q_centroid)
      75             :       CALL section_vals_val_get(normalmode_section, "Q_BEAD", &
      76          54 :                                 r_val=normalmode_env%Q_bead)
      77             :       CALL section_vals_val_get(normalmode_section, "MODEFACTOR", &
      78             :                                 explicit=explicit_modefactor, &
      79          54 :                                 r_val=normalmode_env%modefactor)
      80             :       CALL section_vals_val_get(normalmode_section, "GAMMA", &
      81             :                                 r_val=gamma_parameter, &
      82          54 :                                 explicit=explicit_gamma)
      83             : 
      84          54 :       IF (explicit_modefactor .AND. explicit_gamma) THEN
      85           0 :          CPABORT("Both GAMMA and MODEFACTOR have been declared. Please use only one.")
      86             :       END IF
      87          54 :       IF (explicit_gamma) THEN
      88           4 :          normalmode_env%modefactor = 1.0_dp/gamma_parameter**2
      89             :       END IF
      90             : 
      91          54 :       IF (propagator == propagator_cmd) THEN
      92           4 :          IF (.NOT. explicit_gamma) THEN
      93           0 :             CPABORT("GAMMA needs to be specified with CMD PROPAGATOR")
      94             :          END IF
      95           4 :          IF (gamma_parameter <= 1.0_dp) THEN
      96           0 :             CPWARN("GAMMA should be larger than 1.0 for CMD PROPAGATOR")
      97             :          END IF
      98             :       END IF
      99             : 
     100          54 :       IF (normalmode_env%Q_centroid < 0.0_dp) THEN
     101          54 :          normalmode_env%Q_centroid = -normalmode_env%Q_centroid/(kT*p)
     102             :       END IF
     103          54 :       IF (normalmode_env%Q_bead < 0.0_dp) THEN
     104          54 :          normalmode_env%Q_bead = -normalmode_env%Q_bead/(kT*p)
     105             :       END IF
     106             : 
     107             :       !Use different normal mode transformations depending on the propagator
     108          54 :       IF (propagator == propagator_pimd .OR. propagator == propagator_cmd) THEN
     109             : 
     110          38 :          IF (propagator == propagator_pimd) THEN
     111          34 :             normalmode_env%harm = p*kT*kT/normalmode_env%modefactor
     112           4 :          ELSE IF (propagator == propagator_cmd) THEN
     113           4 :             normalmode_env%harm = p*kT*kT*gamma_parameter*gamma_parameter
     114           4 :             normalmode_env%modefactor = 1.0_dp/(gamma_parameter*gamma_parameter)
     115             :          END IF
     116             : 
     117             :          ! set up the transformation matrices
     118         214 :          DO i = 1, p
     119         176 :             normalmode_env%lambda(i) = 2.0_dp*(1.0_dp - COS(pi*(i/2)*2.0_dp/p))
     120        1206 :             DO j = 1, p
     121         992 :                k = ((i/2)*(j - 1))/p
     122         992 :                k = (i/2)*(j - 1) - k*p
     123         992 :                li = 2*(i - 2*(i/2))*p - p
     124        1168 :                normalmode_env%u2x(j, i) = SQRT(2.0_dp/p)*SIN(twopi*(k + 0.125_dp*li)/p)
     125             :             END DO
     126             :          END DO
     127          38 :          normalmode_env%lambda(1) = 1.0_dp/(p*normalmode_env%modefactor)
     128         214 :          DO i = 1, p
     129        1206 :             DO j = 1, p
     130             :                normalmode_env%x2u(i, j) = SQRT(normalmode_env%lambda(i)* &
     131             :                                                normalmode_env%modefactor)* &
     132        1168 :                                           normalmode_env%u2x(j, i)
     133             :             END DO
     134             :          END DO
     135         214 :          DO i = 1, p
     136        1206 :             DO j = 1, p
     137             :                normalmode_env%u2x(i, j) = normalmode_env%u2x(i, j)/ &
     138             :                                           SQRT(normalmode_env%lambda(j)* &
     139        1168 :                                                normalmode_env%modefactor)
     140             :             END DO
     141             :          END DO
     142         214 :          normalmode_env%lambda(:) = normalmode_env%harm
     143             : 
     144          16 :       ELSE IF (propagator == propagator_rpmd) THEN
     145          16 :          normalmode_env%harm = kT/normalmode_env%modefactor
     146          16 :          sqrt2p = SQRT(2.0_dp/REAL(p, dp))
     147          16 :          twopip = twopi/REAL(p, dp)
     148          16 :          pip = pi/REAL(p, dp)
     149          16 :          invsqrtp = 1.0_dp/SQRT(REAL(p, dp))
     150        2680 :          normalmode_env%x2u(:, :) = 0.0_dp
     151         152 :          normalmode_env%x2u(1, :) = invsqrtp
     152         152 :          DO j = 1, p
     153        1400 :             DO i = 2, p/2 + 1
     154        1400 :                normalmode_env%x2u(i, j) = sqrt2p*COS(twopip*(i - 1)*(j - 1))
     155             :             END DO
     156        1280 :             DO i = p/2 + 2, p
     157        1264 :                normalmode_env%x2u(i, j) = sqrt2p*SIN(twopip*(i - 1)*(j - 1))
     158             :             END DO
     159             :          END DO
     160          16 :          IF (MOD(p, 2) == 0) THEN
     161          16 :             DO i = 1, p - 1, 2
     162          68 :                normalmode_env%x2u(p/2 + 1, i) = invsqrtp
     163          68 :                normalmode_env%x2u(p/2 + 1, i + 1) = -1.0_dp*invsqrtp
     164             :             END DO
     165             :          END IF
     166             : 
     167        5344 :          normalmode_env%u2x = TRANSPOSE(normalmode_env%x2u)
     168             : 
     169             :          ! Setting up propagator frequencies for rpmd
     170          16 :          normalmode_env%lambda(1) = 0.0_dp
     171         136 :          DO i = 2, p
     172         120 :             normalmode_env%lambda(i) = 2.0_dp*normalmode_env%harm*SIN((i - 1)*pip)
     173         136 :             normalmode_env%lambda(i) = normalmode_env%lambda(i)*normalmode_env%lambda(i)
     174             :          END DO
     175          16 :          normalmode_env%harm = kT*kT
     176             :       ELSE
     177           0 :          CPABORT("UNKNOWN PROPAGATOR FOR PINT SELECTED")
     178             :       END IF
     179             : 
     180          54 :    END SUBROUTINE normalmode_env_create
     181             : 
     182             : ! ***************************************************************************
     183             : !> \brief releases the normalmode environment
     184             : !> \param normalmode_env the normalmode_env to release
     185             : !> \author Harald Forbert
     186             : ! **************************************************************************************************
     187          54 :    PURE SUBROUTINE normalmode_release(normalmode_env)
     188             : 
     189             :       TYPE(normalmode_env_type), INTENT(INOUT)           :: normalmode_env
     190             : 
     191          54 :       DEALLOCATE (normalmode_env%x2u)
     192          54 :       DEALLOCATE (normalmode_env%u2x)
     193          54 :       DEALLOCATE (normalmode_env%lambda)
     194             : 
     195          54 :    END SUBROUTINE normalmode_release
     196             : 
     197             : ! ***************************************************************************
     198             : !> \brief initializes the masses and fictitious masses compatible with the
     199             : !>      normal mode information
     200             : !> \param normalmode_env the definition of the normal mode transformation
     201             : !> \param mass *input* the masses of the particles
     202             : !> \param mass_beads masses of the beads
     203             : !> \param mass_fict the fictitious masses
     204             : !> \param Q masses of the nose thermostats
     205             : !> \author Harald Forbert
     206             : ! **************************************************************************************************
     207          54 :    PURE SUBROUTINE normalmode_init_masses(normalmode_env, mass, mass_beads, mass_fict, &
     208          54 :                                           Q)
     209             : 
     210             :       TYPE(normalmode_env_type), INTENT(IN)              :: normalmode_env
     211             :       REAL(kind=dp), DIMENSION(:), INTENT(in)            :: mass
     212             :       REAL(kind=dp), DIMENSION(:, :), INTENT(out), &
     213             :          OPTIONAL                                        :: mass_beads, mass_fict
     214             :       REAL(kind=dp), DIMENSION(:), INTENT(out), OPTIONAL :: Q
     215             : 
     216             :       INTEGER                                            :: iat, ib
     217             : 
     218          54 :       IF (PRESENT(Q)) THEN
     219         366 :          Q = normalmode_env%Q_bead
     220          54 :          Q(1) = normalmode_env%Q_centroid
     221             :       END IF
     222          54 :       IF (PRESENT(mass_beads) .OR. PRESENT(mass_fict)) THEN
     223          54 :          IF (PRESENT(mass_beads)) THEN
     224       64944 :             DO iat = 1, SIZE(mass)
     225       64890 :                mass_beads(1, iat) = 0.0_dp
     226      297342 :                DO ib = 2, normalmode_env%p
     227      297288 :                   mass_beads(ib, iat) = mass(iat)
     228             :                END DO
     229             :             END DO
     230             :          END IF
     231          54 :          IF (PRESENT(mass_fict)) THEN
     232       64944 :             DO iat = 1, SIZE(mass)
     233      362232 :                DO ib = 1, normalmode_env%p
     234      362178 :                   mass_fict(ib, iat) = mass(iat)
     235             :                END DO
     236             :             END DO
     237             :          END IF
     238             :       END IF
     239             : 
     240          54 :    END SUBROUTINE normalmode_init_masses
     241             : 
     242             : ! ***************************************************************************
     243             : !> \brief Transforms from the x into the u variables using a normal mode
     244             : !>      transformation for the positions
     245             : !> \param normalmode_env the environment for the normal mode transformation
     246             : !> \param ux will contain the u variable
     247             : !> \param x the positions to transform
     248             : !> \author Harald Forbert
     249             : ! **************************************************************************************************
     250          91 :    SUBROUTINE normalmode_x2u(normalmode_env, ux, x)
     251             :       TYPE(normalmode_env_type), INTENT(INOUT)           :: normalmode_env
     252             :       REAL(kind=dp), DIMENSION(:, :), INTENT(out)        :: ux
     253             :       REAL(kind=dp), DIMENSION(:, :), INTENT(in)         :: x
     254             : 
     255             :       CALL DGEMM('N', 'N', normalmode_env%p, SIZE(x, 2), normalmode_env%p, 1.0_dp, &
     256             :                  normalmode_env%x2u(1, 1), SIZE(normalmode_env%x2u, 1), x(1, 1), SIZE(x, 1), &
     257          91 :                  0.0_dp, ux, SIZE(ux, 1))
     258          91 :    END SUBROUTINE normalmode_x2u
     259             : 
     260             : ! ***************************************************************************
     261             : !> \brief transform from the u variable to the x (back normal mode
     262             : !>      transformation for the positions)
     263             : !> \param normalmode_env the environment for the normal mode transformation
     264             : !> \param ux the u variable (positions to be backtransformed)
     265             : !> \param x will contain the positions
     266             : !> \author Harald Forbert
     267             : ! **************************************************************************************************
     268        2265 :    SUBROUTINE normalmode_u2x(normalmode_env, ux, x)
     269             :       TYPE(normalmode_env_type), INTENT(INOUT)           :: normalmode_env
     270             :       REAL(kind=dp), DIMENSION(:, :), INTENT(in)         :: ux
     271             :       REAL(kind=dp), DIMENSION(:, :), INTENT(out)        :: x
     272             : 
     273             :       CALL DGEMM('N', 'N', normalmode_env%p, SIZE(ux, 2), normalmode_env%p, 1.0_dp, &
     274             :                  normalmode_env%u2x(1, 1), SIZE(normalmode_env%u2x, 1), ux(1, 1), SIZE(ux, 1), &
     275        2265 :                  0.0_dp, x, SIZE(x, 1))
     276        2265 :    END SUBROUTINE normalmode_u2x
     277             : 
     278             : ! ***************************************************************************
     279             : !> \brief normalmode transformation for the forces
     280             : !> \param normalmode_env the environment for the normal mode transformation
     281             : !> \param uf will contain the forces for the transformed variables afterwards
     282             : !> \param f the forces to transform
     283             : !> \author Harald Forbert
     284             : ! **************************************************************************************************
     285         496 :    SUBROUTINE normalmode_f2uf(normalmode_env, uf, f)
     286             :       TYPE(normalmode_env_type), INTENT(INOUT)           :: normalmode_env
     287             :       REAL(kind=dp), DIMENSION(:, :), INTENT(out)        :: uf
     288             :       REAL(kind=dp), DIMENSION(:, :), INTENT(in)         :: f
     289             : 
     290             :       CALL DGEMM('T', 'N', normalmode_env%p, SIZE(f, 2), normalmode_env%p, 1.0_dp, &
     291             :                  normalmode_env%u2x(1, 1), SIZE(normalmode_env%u2x, 1), f(1, 1), SIZE(f, 1), &
     292         496 :                  0.0_dp, uf, SIZE(uf, 1))
     293         496 :    END SUBROUTINE normalmode_f2uf
     294             : 
     295             : ! ***************************************************************************
     296             : !> \brief calculates the harmonic force in the normal mode basis
     297             : !> \param normalmode_env the normal mode environment
     298             : !> \param mass_beads the masses of the beads
     299             : !> \param ux the positions of the beads in the staging basis
     300             : !> \param uf_h the harmonic forces (not accelerations)
     301             : !> \param e_h ...
     302             : !> \author Harald Forbert
     303             : ! **************************************************************************************************
     304        1226 :    PURE SUBROUTINE normalmode_calc_uf_h(normalmode_env, mass_beads, ux, uf_h, e_h)
     305             :       TYPE(normalmode_env_type), INTENT(IN)              :: normalmode_env
     306             :       REAL(kind=dp), DIMENSION(:, :), POINTER            :: mass_beads, ux, uf_h
     307             :       REAL(KIND=dp), INTENT(OUT)                         :: e_h
     308             : 
     309             :       INTEGER                                            :: ibead, idim
     310             :       REAL(kind=dp)                                      :: f
     311             : 
     312        1226 :       e_h = 0.0_dp
     313      446582 :       DO idim = 1, SIZE(mass_beads, 2)
     314             : 
     315             :          ! starting at 2 since the centroid is at 1 and it's mass_beads
     316             :          ! SHOULD be zero anyways:
     317             : 
     318      445356 :          uf_h(1, idim) = 0.0_dp
     319     2114426 :          DO ibead = 2, normalmode_env%p
     320     1667844 :             f = -mass_beads(ibead, idim)*normalmode_env%lambda(ibead)*ux(ibead, idim)
     321     1667844 :             uf_h(ibead, idim) = f
     322             :             ! - to cancel the - in the force f.
     323     2113200 :             e_h = e_h - 0.5_dp*ux(ibead, idim)*f
     324             :          END DO
     325             : 
     326             :       END DO
     327        1226 :    END SUBROUTINE normalmode_calc_uf_h
     328             : 
     329             : END MODULE pint_normalmode