LCOV - code coverage report
Current view: top level - src/motion - pint_piglet.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:5911e77) Lines: 208 263 79.1 %
Date: 2024-12-12 07:20:04 Functions: 9 9 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief  Methods to apply the piglet thermostat to PI runs.
      10             : !> \author Felix Uhl
      11             : !> \par History
      12             : !>      10.2014 created [Felix Uhl]
      13             : ! **************************************************************************************************
      14             : MODULE pint_piglet
      15             :    USE cp_files,                        ONLY: close_file,&
      16             :                                               open_file
      17             :    USE cp_units,                        ONLY: cp_unit_from_cp2k,&
      18             :                                               cp_unit_to_cp2k
      19             :    USE gle_system_dynamics,             ONLY: gle_cholesky_stab,&
      20             :                                               gle_matrix_exp
      21             :    USE input_constants,                 ONLY: matrix_init_cholesky,&
      22             :                                               matrix_init_diagonal,&
      23             :                                               propagator_rpmd
      24             :    USE input_section_types,             ONLY: section_vals_get,&
      25             :                                               section_vals_get_subs_vals,&
      26             :                                               section_vals_type,&
      27             :                                               section_vals_val_get
      28             :    USE kinds,                           ONLY: default_path_length,&
      29             :                                               default_string_length,&
      30             :                                               dp
      31             :    USE message_passing,                 ONLY: mp_para_env_type
      32             :    USE parallel_rng_types,              ONLY: GAUSSIAN,&
      33             :                                               rng_record_length,&
      34             :                                               rng_stream_type,&
      35             :                                               rng_stream_type_from_record
      36             :    USE pint_io,                         ONLY: pint_write_line
      37             :    USE pint_types,                      ONLY: piglet_therm_type,&
      38             :                                               pint_env_type
      39             : #include "../base/base_uses.f90"
      40             : 
      41             :    IMPLICIT NONE
      42             : 
      43             :    PRIVATE
      44             : 
      45             :    PUBLIC :: pint_piglet_create, &
      46             :              pint_piglet_init, &
      47             :              pint_piglet_step, &
      48             :              pint_piglet_release, &
      49             :              pint_calc_piglet_energy
      50             : 
      51             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pint_piglet'
      52             : 
      53             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      54             : !                                           !
      55             : !       _._ _..._ .-',     _.._   _         !
      56             : !      '-. `     '  /-._.-'    ',/ )        !
      57             : !         )         \            '.         !
      58             : !        / _    _    |             \        !
      59             : !       |  o    o    /              |       !
      60             : !       \   .-.                     ;       !
      61             : !        '-('' ).-'       ,'       ;        !
      62             : !           '-;           |      .'         !
      63             : !              \           \    /           !
      64             : !              | 7  .__  _.-\   \           !
      65             : !              | |  |  ``/  /`  /           !
      66             : !             /,_|  |   /,_/   /            !
      67             : !                /,_/      '`-'             !
      68             : !                                           !
      69             : !    (    (             (                   !
      70             : !    )\ ) )\ )  (       )\ )        *   )   !
      71             : !   (()/((()/(  )\ )   (()/(  (   ` )  /(   !
      72             : !    /(_))/(_))(()/(    /(_)) )\   ( )(_))  !
      73             : !   (_)) (_))   /(_))_ (_))  ((_) (_(_())   !
      74             : !   | _ \|_ _| (_)) __|| |   | __||_   _|   !
      75             : !   |  _/ | |    | (_ || |__ | _|   | |     !
      76             : !   |_|  |___|    \___||____||___|  |_|     !
      77             : !                                           !
      78             : !        Make Quantum Mechanics Hot         !
      79             : !                                           !
      80             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      81             : 
      82             : CONTAINS
      83             : 
      84             : ! ***************************************************************************
      85             : !> \brief creates the data structure for a piglet thermostating in PI runs
      86             : !> \param piglet_therm ...
      87             : !> \param pint_env ...
      88             : !> \param section ...
      89             : !> \author Felix Uhl
      90             : ! **************************************************************************************************
      91          54 :    SUBROUTINE pint_piglet_create(piglet_therm, pint_env, section)
      92             :       TYPE(piglet_therm_type), INTENT(OUT)               :: piglet_therm
      93             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
      94             :       TYPE(section_vals_type), POINTER                   :: section
      95             : 
      96             :       INTEGER                                            :: ndim, p
      97             : 
      98           2 :       CALL section_vals_val_get(section, "NEXTRA_DOF", i_val=piglet_therm%nsp1)
      99             :       !add real degree of freedom to ns to reach nsp1
     100           2 :       piglet_therm%nsp1 = piglet_therm%nsp1 + 1
     101           2 :       p = pint_env%p
     102           2 :       piglet_therm%p = pint_env%p
     103           2 :       ndim = pint_env%ndim
     104           2 :       piglet_therm%ndim = pint_env%ndim
     105           2 :       NULLIFY (piglet_therm%a_mat)
     106          10 :       ALLOCATE (piglet_therm%a_mat(piglet_therm%nsp1, piglet_therm%nsp1, p))
     107        2186 :       piglet_therm%a_mat(:, :, :) = 0.0_dp
     108           2 :       NULLIFY (piglet_therm%c_mat)
     109          10 :       ALLOCATE (piglet_therm%c_mat(piglet_therm%nsp1, piglet_therm%nsp1, p))
     110        2186 :       piglet_therm%c_mat(:, :, :) = 0.0_dp
     111           2 :       NULLIFY (piglet_therm%gle_t)
     112          10 :       ALLOCATE (piglet_therm%gle_t(piglet_therm%nsp1, piglet_therm%nsp1, p))
     113        2186 :       piglet_therm%gle_t(:, :, :) = 0.0_dp
     114           2 :       NULLIFY (piglet_therm%gle_s)
     115          10 :       ALLOCATE (piglet_therm%gle_s(piglet_therm%nsp1, piglet_therm%nsp1, p))
     116        2186 :       piglet_therm%gle_s(:, :, :) = 0.0_dp
     117           2 :       NULLIFY (piglet_therm%smalls)
     118           8 :       ALLOCATE (piglet_therm%smalls(piglet_therm%nsp1, ndim*p))
     119     1105922 :       piglet_therm%smalls(:, :) = 0.0_dp
     120           2 :       NULLIFY (piglet_therm%temp1)
     121           8 :       ALLOCATE (piglet_therm%temp1(piglet_therm%nsp1, ndim))
     122       92162 :       piglet_therm%temp1(:, :) = 0.0_dp
     123           2 :       NULLIFY (piglet_therm%temp2)
     124           8 :       ALLOCATE (piglet_therm%temp2(piglet_therm%nsp1, ndim))
     125       92162 :       piglet_therm%temp2(:, :) = 0.0_dp
     126           2 :       NULLIFY (piglet_therm%sqrtmass)
     127           8 :       ALLOCATE (piglet_therm%sqrtmass(piglet_therm%p, piglet_therm%ndim))
     128      119810 :       piglet_therm%sqrtmass(:, :) = 0.0_dp
     129             : 
     130           2 :    END SUBROUTINE pint_piglet_create
     131             : 
     132             : ! ***************************************************************************
     133             : !> \brief initializes the data for a piglet run
     134             : !> \param piglet_therm ...
     135             : !> \param pint_env ...
     136             : !> \param section ...
     137             : !> \param dt ...
     138             : !> \param para_env ...
     139             : !> \author Felix Uhl
     140             : ! **************************************************************************************************
     141           2 :    SUBROUTINE pint_piglet_init(piglet_therm, pint_env, section, dt, para_env)
     142             :       TYPE(piglet_therm_type), INTENT(INOUT)             :: piglet_therm
     143             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
     144             :       TYPE(section_vals_type), POINTER                   :: section
     145             :       REAL(KIND=dp), INTENT(IN)                          :: dt
     146             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     147             : 
     148             :       CHARACTER(len=20)                                  :: default_format, read_unit
     149             :       CHARACTER(len=default_path_length)                 :: matrices_file_name
     150             :       CHARACTER(len=default_string_length)               :: msg, temp_input
     151             :       CHARACTER(LEN=rng_record_length)                   :: rng_record
     152             :       INTEGER                                            :: cbrac, i, ibead, idim, imode, &
     153             :                                                             input_unit, isp, j, matrix_init, ns, &
     154             :                                                             obrac, p, read_err
     155             :       LOGICAL                                            :: explicit
     156             :       REAL(KIND=dp), DIMENSION(3, 2)                     :: initial_seed
     157           2 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: smallstmp
     158             :       REAL(KIND=dp), &
     159           4 :          DIMENSION(piglet_therm%nsp1, piglet_therm%nsp1) :: Mtmp, tmpmatrix
     160             :       TYPE(section_vals_type), POINTER                   :: rng_section, smalls_section
     161             : 
     162           2 :       p = piglet_therm%p
     163           2 :       pint_env%e_piglet = 0.0_dp
     164           2 :       pint_env%piglet_therm%thermostat_energy = 0.0_dp
     165           2 :       CALL section_vals_val_get(section, "THERMOSTAT_ENERGY", r_val=piglet_therm%thermostat_energy)
     166           2 :       CALL section_vals_val_get(section, "SMATRIX_INIT", i_val=matrix_init)
     167             :       !Read the matices
     168             : 
     169           2 :       IF (pint_env%propagator%prop_kind /= propagator_rpmd) THEN
     170           0 :          CPABORT("PIGLET is designed to work with the RPMD propagator")
     171             :       END IF
     172             :       ! Select algorithm for S-matrix initialization
     173           2 :       IF (matrix_init == matrix_init_cholesky) THEN
     174           0 :          IF (para_env%is_source()) THEN
     175           0 :             CALL pint_write_line("PIGLET| Initalizing S-matrices using cholesky decomposition.")
     176             :          END IF
     177           2 :       ELSE IF (matrix_init == matrix_init_diagonal) THEN
     178           2 :          IF (para_env%is_source()) THEN
     179           1 :             CALL pint_write_line("PIGLET| Initalizing S-matrices using full diagonalization.")
     180             :          END IF
     181             :       ELSE
     182           0 :          CPWARN("No PIGLET init algorithm selected. Selecting cholesky decomposition")
     183           0 :          matrix_init = matrix_init_cholesky
     184             :       END IF
     185             : 
     186           2 :       IF (para_env%is_source()) THEN
     187             :          ! Read input matrices
     188           1 :          WRITE (default_format, '(A,I10,A)') "(A", default_string_length, ")"
     189             :          CALL section_vals_val_get(section, "MATRICES_FILE_NAME", &
     190           1 :                                    c_val=matrices_file_name)
     191           1 :          CALL pint_write_line("PIGLET| Reading PIGLET matrices from file: ")
     192           1 :          CALL pint_write_line("PIGLET|    "//TRIM(matrices_file_name))
     193             :          CALL open_file(file_name=TRIM(matrices_file_name), &
     194           1 :                         file_action="READ", file_status="OLD", unit_number=input_unit)
     195           1 :          read_err = 0
     196           1 :          msg = ""
     197          20 :          DO WHILE (read_err == 0)
     198          20 :             READ (input_unit, default_format, iostat=read_err) temp_input
     199          20 :             IF (read_err /= 0) EXIT
     200             :             !Parse comment section
     201          20 :             IF (INDEX(temp_input, "#") /= 0) THEN
     202          16 :                IF (INDEX(temp_input, "T=") /= 0) THEN
     203             :                   CALL check_temperature(line=temp_input, &
     204             :                                          propagator=pint_env%propagator%prop_kind, &
     205           1 :                                          targettemp=pint_env%kT*pint_env%propagator%temp_sim2phys)
     206          15 :                ELSE IF (INDEX(temp_input, "A MATRIX") /= 0) THEN
     207           1 :                   obrac = INDEX(temp_input, "(") + 1
     208           1 :                   cbrac = INDEX(temp_input, ")") - 1
     209           1 :                   read_unit = temp_input(obrac:cbrac)
     210          13 :                   DO imode = 1, p
     211          12 :                      READ (input_unit, default_format) temp_input
     212         121 :                      DO i = 1, piglet_therm%nsp1
     213             :                         READ (input_unit, *, iostat=read_err) &
     214        1080 :                            (piglet_therm%a_mat(i, j, imode), j=1, piglet_therm%nsp1)
     215         120 :                         IF (read_err /= 0) THEN
     216           0 :                            WRITE (UNIT=msg, FMT=*) "Invalid PIGLET A-matrix Nr.", i - 1
     217           0 :                            CPABORT(msg)
     218           0 :                            EXIT
     219             :                         END IF
     220             :                      END DO
     221             :                   END DO
     222             :                   !convert to cp2k units
     223           1 :                   IF (read_err == 0) THEN
     224             :                      CALL a_mat_to_cp2k(piglet_therm%a_mat, p, &
     225           1 :                                         piglet_therm%nsp1, read_unit, msg)
     226             :                   END IF
     227          14 :                ELSE IF (INDEX(temp_input, "C MATRIX") /= 0) THEN
     228           1 :                   obrac = INDEX(temp_input, "(") + 1
     229           1 :                   cbrac = INDEX(temp_input, ")") - 1
     230           1 :                   read_unit = temp_input(obrac:cbrac)
     231          13 :                   DO imode = 1, p
     232          12 :                      READ (input_unit, default_format) temp_input
     233         121 :                      DO i = 1, piglet_therm%nsp1
     234             :                         READ (input_unit, *, iostat=read_err) &
     235        1080 :                            (piglet_therm%c_mat(i, j, imode), j=1, piglet_therm%nsp1)
     236         120 :                         IF (read_err /= 0) THEN
     237           0 :                            WRITE (UNIT=msg, FMT=*) "Invalid PIGLET C-matrix Nr.", i - 1
     238           0 :                            CPABORT(msg)
     239           0 :                            EXIT
     240             :                         END IF
     241             :                      END DO
     242             :                   END DO
     243           1 :                   IF (read_err == 0) THEN
     244             :                      CALL c_mat_to_cp2k(piglet_therm%c_mat, p, &
     245           1 :                                         piglet_therm%nsp1, read_unit, msg)
     246             :                   END IF
     247             :                END IF
     248             :             END IF
     249             :          END DO
     250           1 :          CALL close_file(unit_number=input_unit)
     251             :       END IF
     252             :       ! communicate A and C matrix to other nodes
     253             :       CALL para_env%bcast(piglet_therm%a_mat, &
     254        4370 :                           para_env%source)
     255             :       CALL para_env%bcast(piglet_therm%c_mat, &
     256        4370 :                           para_env%source)
     257             : 
     258             :       !prepare Random number generator
     259           2 :       NULLIFY (rng_section)
     260             :       rng_section => section_vals_get_subs_vals(section, &
     261           2 :                                                 subsection_name="RNG_INIT")
     262           2 :       CALL section_vals_get(rng_section, explicit=explicit)
     263           2 :       IF (explicit) THEN
     264             :          CALL section_vals_val_get(rng_section, "_DEFAULT_KEYWORD_", &
     265           0 :                                    i_rep_val=1, c_val=rng_record)
     266           0 :          piglet_therm%gaussian_rng_stream = rng_stream_type_from_record(rng_record)
     267             :       ELSE
     268          18 :          initial_seed(:, :) = REAL(pint_env%thermostat_rng_seed, dp)
     269             :          piglet_therm%gaussian_rng_stream = rng_stream_type( &
     270             :                                             name="piglet_rng_gaussian", distribution_type=GAUSSIAN, &
     271             :                                             extended_precision=.TRUE., &
     272           2 :                                             seed=initial_seed)
     273             :       END IF
     274             : 
     275             :       !Compute the T and S matrices on every mpi process
     276          26 :       DO i = 1, p
     277             :          ! T = EXP(-A_mat*dt) = EXP(-A_mat*dt/2)
     278             :          ! Values for j and k = 15 are way to high, but to be sure.
     279             :          ! (its only executed once anyway)
     280             :          CALL gle_matrix_exp(M=(-0.5_dp*dt)*piglet_therm%a_mat(:, :, i), & !dt scaled A matrix
     281             :                              n=piglet_therm%nsp1, & !size of matrix
     282             :                              j=15, & !truncation for taylor expansion
     283             :                              k=15, & !scaling parameter for faster convergence of taylor expansion
     284        2184 :                              EM=piglet_therm%gle_t(:, :, i)) !output T matrices
     285             :          ! S*TRANSPOSE(S) = C-T*C*TRANSPOSE(T)
     286             :          ! T*C:
     287             :          CALL DGEMM('N', & ! T-matrix is not transposed
     288             :                     'N', & ! C-matrix is not transposed
     289             :                     piglet_therm%nsp1, & ! number of rows of T-matrix
     290             :                     piglet_therm%nsp1, & ! number of columns of C-matrix
     291             :                     piglet_therm%nsp1, & ! number of columns of T-matrix and number of rows of C-matrix
     292             :                     1.0D0, & ! scaling factor alpha
     293             :                     piglet_therm%gle_t(:, :, i), & ! T-matrix
     294             :                     piglet_therm%nsp1, & ! leading dimension of T-matrix as declared
     295             :                     piglet_therm%c_mat(:, :, i), & ! C-matrix
     296             :                     piglet_therm%nsp1, & ! leading dimension of C-matrix as declared
     297             :                     0.0D0, & ! scaling of tmpmatrix as additive
     298             :                     tmpmatrix, & ! result matrix: tmpmatrix
     299          24 :                     piglet_therm%nsp1) ! leading dimension of tmpmatrix
     300             :          ! T*C*TRANSPOSE(T):
     301             :          CALL DGEMM('N', & ! tmpmatrix is not transposed
     302             :                     'T', & ! T-matrix is transposed
     303             :                     piglet_therm%nsp1, & ! number of rows of tmpmatrix
     304             :                     piglet_therm%nsp1, & ! number of columns of T-matrix
     305             :                     piglet_therm%nsp1, & ! number of columns of tmpmatrix and number of rows of T-matrix
     306             :                     1.0D0, & ! scaling factor alpha
     307             :                     tmpmatrix, & ! tmpmatrix
     308             :                     piglet_therm%nsp1, & ! leading dimension of tmpmatrix as declared
     309             :                     piglet_therm%gle_t(:, :, i), & ! T-matrix
     310             :                     piglet_therm%nsp1, & ! leading dimension of T-matrix as declared
     311             :                     0.0D0, & ! scaling of Mtmp as additive
     312             :                     Mtmp, & ! result matrix: Mtmp
     313          24 :                     piglet_therm%nsp1) ! leading dimension of Mtmp
     314             :          ! C - T*C*TRANSPOSE(T):
     315        2184 :          Mtmp(:, :) = piglet_therm%c_mat(:, :, i) - Mtmp(:, :)
     316             : 
     317          26 :          IF (matrix_init == matrix_init_cholesky) THEN
     318             :             ! Get S by cholesky decomposition of Mtmp
     319             :             CALL gle_cholesky_stab(Mtmp, & ! Matrix to decompose
     320             :                                    piglet_therm%gle_s(:, :, i), & ! result
     321           0 :                                    piglet_therm%nsp1) ! Size of the matrix
     322          24 :          ELSE IF (matrix_init == matrix_init_diagonal) THEN
     323             :             ! Get S by full diagonalization of MTmp matrix
     324             :             CALL sqrt_pos_def_mat(piglet_therm%nsp1, & ! Size of the matrix
     325             :                                   Mtmp, & ! matrix to decompose
     326          24 :                                   piglet_therm%gle_s(:, :, i)) ! result
     327             :          END IF
     328             : 
     329             :       END DO
     330             : 
     331             :       ! Initialize extra degrees of freedom for Markovian Dynamics
     332             :       ! as a cholesky decomposition of C-matrix multiplied by a random number vector
     333             :       ! Or from restart
     334     1105922 :       piglet_therm%smalls = 0.0_dp
     335           2 :       NULLIFY (smalls_section)
     336           2 :       smalls_section => section_vals_get_subs_vals(section, subsection_name="EXTRA_DOF")
     337           2 :       CALL section_vals_get(smalls_section, explicit=explicit)
     338           2 :       IF (explicit) THEN
     339           0 :          NULLIFY (smallstmp)
     340             :          CALL section_vals_val_get(smalls_section, "_DEFAULT_KEYWORD_", &
     341           0 :                                    n_rep_val=ns)
     342             :          CALL section_vals_val_get(smalls_section, "_DEFAULT_KEYWORD_", &
     343           0 :                                    r_vals=smallstmp)
     344           0 :          i = 1
     345           0 :          DO isp = 2, piglet_therm%nsp1
     346           0 :             DO ibead = 1, piglet_therm%p*piglet_therm%ndim
     347           0 :                piglet_therm%smalls(isp, ibead) = smallstmp(i)
     348           0 :                i = i + 1
     349             :             END DO
     350             :          END DO
     351             :       ELSE
     352          26 :          DO ibead = 1, piglet_therm%p
     353          24 :             IF (matrix_init == matrix_init_cholesky) THEN
     354             :                CALL gle_cholesky_stab(piglet_therm%c_mat(:, :, ibead), & ! Matrix to decompose
     355             :                                       Mtmp, & ! Result
     356           0 :                                       piglet_therm%nsp1) ! Size of Matrix
     357          24 :             ELSE IF (matrix_init == matrix_init_diagonal) THEN
     358             :                ! Get S by full diagonalization of c_mat matrix
     359             :                CALL sqrt_pos_def_mat(piglet_therm%nsp1, & ! Size of the matrix
     360             :                                      piglet_therm%c_mat(:, :, ibead), & ! matrix to decompose
     361          24 :                                      Mtmp) ! result
     362             :             END IF
     363             :             ! Fill a vector with random numbers
     364      110616 :             DO idim = 1, piglet_therm%ndim
     365     1105944 :                DO j = 1, piglet_therm%nsp1
     366     1105920 :                   piglet_therm%temp2(j, idim) = piglet_therm%gaussian_rng_stream%next()
     367             :                   !piglet_therm%temp2(j,idim) = 1.0_dp
     368             :                END DO
     369             :             END DO
     370             :             CALL DGEMM("N", & ! Matrix Mtmp is not transposed
     371             :                        "N", & ! Matrix temp2 is not transposed
     372             :                        piglet_therm%nsp1, & ! Number of rows of matrix Mtmp
     373             :                        piglet_therm%ndim, & ! Number of columns of matrix temp2
     374             :                        piglet_therm%nsp1, & ! Number of columns of matrix Mtmp
     375             :                        1.0_dp, & !scaling of Mtmp
     376             :                        Mtmp(1, 1), & ! Matrix Mtmp
     377             :                        piglet_therm%nsp1, & ! leading dimension of Mtmp
     378             :                        piglet_therm%temp2, & ! temp2 matrix
     379             :                        piglet_therm%nsp1, & ! leading dimension of temp2
     380             :                        0.0_dp, & ! scaling of added matrix smalls
     381             :                        piglet_therm%temp1, & ! result matrix
     382          24 :                        piglet_therm%nsp1) ! leading dimension of result matrix
     383             : 
     384      110618 :             DO idim = 1, piglet_therm%ndim
     385      110592 :                j = (idim - 1)*piglet_therm%p + ibead
     386     1105944 :                DO i = 1, piglet_therm%nsp1
     387     1105920 :                   piglet_therm%smalls(i, j) = piglet_therm%temp1(i, idim)
     388             :                END DO
     389             :             END DO
     390             :          END DO
     391             :       END IF
     392             : 
     393             :       !Fill the array for the sqrt of the masses
     394        9218 :       DO idim = 1, piglet_therm%ndim
     395      119810 :          DO ibead = 1, piglet_therm%p
     396      119808 :             piglet_therm%sqrtmass(ibead, idim) = SQRT(pint_env%mass_fict(ibead, idim))
     397             :          END DO
     398             :       END DO
     399             : 
     400           2 :    END SUBROUTINE pint_piglet_init
     401             : 
     402             : ! **************************************************************************************************
     403             : !> \brief ...
     404             : !> \param vold ...
     405             : !> \param vnew ...
     406             : !> \param first_mode ...
     407             : !> \param masses ...
     408             : !> \param piglet_therm ...
     409             : ! **************************************************************************************************
     410          20 :    SUBROUTINE pint_piglet_step(vold, vnew, first_mode, masses, piglet_therm)
     411             :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: vold, vnew
     412             :       INTEGER, INTENT(IN)                                :: first_mode
     413             :       REAL(kind=dp), DIMENSION(:, :), INTENT(IN)         :: masses
     414             :       TYPE(piglet_therm_type), POINTER                   :: piglet_therm
     415             : 
     416             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'pint_piglet_step'
     417             : 
     418             :       INTEGER                                            :: handle, i, ibead, idim, j, ndim, ndimp, &
     419             :                                                             nsp1, p
     420             :       REAL(KIND=dp)                                      :: delta_ekin
     421             : 
     422          20 :       CALL timeset(routineN, handle)
     423          20 :       nsp1 = piglet_therm%nsp1
     424          20 :       ndim = piglet_therm%ndim
     425          20 :       p = piglet_therm%p
     426          20 :       ndimp = ndim*p
     427             :       ! perform the following operation for all 3N*P
     428             :       ! smalls = gle_t*smalls + gle_s*rand_mat
     429             :       ! Copy mass scaled momenta to temp1 matrix
     430             :       ! p/sqrt(m) we use velocities so v*sqrt(m)
     431         260 :       DO ibead = first_mode, p
     432             :          ! Copy mass scaled momenta to temp1 matrix
     433             :          ! p/sqrt(m) we use velocities so v*sqrt(m)
     434     1106160 :          DO idim = 1, ndim
     435     1106160 :             piglet_therm%temp1(1, idim) = vold(ibead, idim)*piglet_therm%sqrtmass(ibead, idim)
     436             :          END DO
     437             :          ! copy the extra degrees of freedom to the temp1 matrix
     438     1106160 :          DO idim = 1, ndim
     439     9953520 :             DO i = 2, nsp1
     440     9953280 :                piglet_therm%temp1(i, idim) = piglet_therm%smalls(i, (ibead - 1)*ndim + idim)
     441             :             END DO
     442             :          END DO
     443             : 
     444             :          !fill temp2 with gaussian random noise
     445        2400 :          DO j = 1, nsp1
     446     9955680 :             DO idim = 1, ndim
     447     9955440 :                piglet_therm%temp2(j, idim) = piglet_therm%gaussian_rng_stream%next()
     448             :                !piglet_therm%temp2(j,idim) = 1.0_dp
     449             :             END DO
     450             :          END DO
     451             : 
     452         240 :          i = (ibead - 1)*piglet_therm%ndim + 1
     453             :          !smalls(:,i) = 1*S*temp2 + 0 * smalls
     454             :          CALL DGEMM("N", & ! S-matrix should not be transposed
     455             :                     "N", & ! tmp2 matrix shoud not be transposed
     456             :                     nsp1, & ! Number of rows of S-Matrix
     457             :                     ndim, & ! Number of columns of temp2 vector
     458             :                     nsp1, & ! Number of Columns of S-Matrix
     459             :                     1.0_dp, & ! Scaling of S-Matrix
     460             :                     piglet_therm%gle_s(:, :, ibead), & ! S-matrix
     461             :                     nsp1, & ! Leading dimension of S-matrix
     462             :                     piglet_therm%temp2, & ! temp2 vector
     463             :                     nsp1, & ! Leading dimension of temp2
     464             :                     0.0_dp, & ! scaling factor of added smalls vector
     465             :                     piglet_therm%smalls(:, i), & ! result vector
     466         240 :                     nsp1) ! Leading dimension of result vector
     467             : 
     468             :          ! Now add the product of T-matrix * old smalls vectors
     469             :          ! smalls (:,i) = 1*T*temp1 + 1*smalls
     470             :          CALL DGEMM("N", & ! T-matrix should not be transposed
     471             :                     "N", & ! temp1 matrix shoud not be transposed
     472             :                     nsp1, & ! Number of rows of T-Matrix
     473             :                     ndim, & ! Number of columns of temp1 vector
     474             :                     nsp1, & ! Number of Columns of T-Matrix
     475             :                     1.0_dp, & ! Scaling of T-Matrix
     476             :                     piglet_therm%gle_t(:, :, ibead), & ! T-matrix
     477             :                     nsp1, & ! Leading dimension of T-matrix
     478             :                     piglet_therm%temp1, & ! temp1 vector
     479             :                     nsp1, & ! Leading dimension of temp1
     480             :                     1.0_dp, & ! scaling factor of added smalls vector
     481             :                     piglet_therm%smalls(:, i), & ! result vector
     482         260 :                     nsp1) ! Leading dimension of result vector
     483             :       END DO
     484             : 
     485             :       ! Copy the mass scales momenta to the outgoing velocities
     486          20 :       delta_ekin = 0.0_dp
     487       92180 :       DO idim = 1, ndim
     488     1198100 :          DO ibead = 1, p
     489     1105920 :             vnew(ibead, idim) = piglet_therm%smalls(1, (ibead - 1)*ndim + idim)/piglet_therm%sqrtmass(ibead, idim)
     490             :             delta_ekin = delta_ekin + masses(ibead, idim)*( &
     491             :                          vnew(ibead, idim)*vnew(ibead, idim) - &
     492     1198080 :                          vold(ibead, idim)*vold(ibead, idim))
     493             :          END DO
     494             :       END DO
     495             : 
     496             :       ! the piglet is such a strong thermostat, that it messes up the "exact" integration. The thermostats energy will rise lineary, because "it will suck up its own mess" (quote from Michele Ceriotti)
     497          20 :       piglet_therm%thermostat_energy = piglet_therm%thermostat_energy - 0.5_dp*delta_ekin
     498             : 
     499          20 :       CALL timestop(handle)
     500             : 
     501          20 :    END SUBROUTINE pint_piglet_step
     502             : 
     503             : ! ***************************************************************************
     504             : !> \brief releases the piglet environment
     505             : !> \param piglet_therm piglet data to be released
     506             : !> \author Felix Uhl
     507             : ! **************************************************************************************************
     508           2 :    SUBROUTINE pint_piglet_release(piglet_therm)
     509             : 
     510             :       TYPE(piglet_therm_type), INTENT(INOUT)             :: piglet_therm
     511             : 
     512           2 :       DEALLOCATE (piglet_therm%a_mat)
     513           2 :       DEALLOCATE (piglet_therm%c_mat)
     514           2 :       DEALLOCATE (piglet_therm%gle_t)
     515           2 :       DEALLOCATE (piglet_therm%gle_s)
     516           2 :       DEALLOCATE (piglet_therm%smalls)
     517           2 :       DEALLOCATE (piglet_therm%temp1)
     518           2 :       DEALLOCATE (piglet_therm%temp2)
     519           2 :       DEALLOCATE (piglet_therm%sqrtmass)
     520             : 
     521           2 :    END SUBROUTINE pint_piglet_release
     522             : 
     523             : ! ***************************************************************************
     524             : !> \brief adjust the unit of A MAT for piglet
     525             : !> \param a_mat ...
     526             : !> \param p ...
     527             : !> \param nsp1 ...
     528             : !> \param myunit ...
     529             : !> \param msg ...
     530             : !> \author Felix Uhl
     531             : ! **************************************************************************************************
     532           1 :    SUBROUTINE a_mat_to_cp2k(a_mat, p, nsp1, myunit, msg)
     533             : 
     534             :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: a_mat
     535             :       INTEGER, INTENT(IN)                                :: p, nsp1
     536             :       CHARACTER(LEN=20), INTENT(IN)                      :: myunit
     537             :       CHARACTER(default_string_length), INTENT(OUT)      :: msg
     538             : 
     539             :       CHARACTER(LEN=20)                                  :: isunit
     540             :       INTEGER                                            :: i, imode, j
     541             : 
     542           1 :       msg = ""
     543           1 :       SELECT CASE (TRIM(myunit))
     544             :       CASE ("femtoseconds^-1")
     545           0 :          isunit = "fs^-1"
     546             :       CASE ("picoseconds^-1")
     547           1 :          isunit = "ps^-1"
     548             :       CASE ("seconds^-1")
     549           0 :          isunit = "s^-1"
     550             :       CASE ("atomic time units^-1")
     551           0 :          RETURN
     552             :       CASE DEFAULT
     553           0 :          msg = "Unknown unit of A-Matrices for PIGLET. Assuming a.u."
     554           0 :          CPWARN(msg)
     555           1 :          RETURN
     556             :       END SELECT
     557             : 
     558          13 :       DO imode = 1, p
     559         121 :          DO j = 1, nsp1
     560        1092 :             DO i = 1, nsp1
     561        1080 :                a_mat(i, j, imode) = cp_unit_to_cp2k(a_mat(i, j, imode), TRIM(isunit))
     562             :             END DO
     563             :          END DO
     564             :       END DO
     565             : 
     566             :    END SUBROUTINE
     567             : ! ***************************************************************************
     568             : !> \brief adjust the unit of C MAT for piglet
     569             : !> \param c_mat ...
     570             : !> \param p ...
     571             : !> \param nsp1 ...
     572             : !> \param myunit ...
     573             : !> \param msg ...
     574             : !> \author Felix Uhl
     575             : ! **************************************************************************************************
     576           1 :    SUBROUTINE c_mat_to_cp2k(c_mat, p, nsp1, myunit, msg)
     577             : 
     578             :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: c_mat
     579             :       INTEGER, INTENT(IN)                                :: p, nsp1
     580             :       CHARACTER(LEN=20), INTENT(IN)                      :: myunit
     581             :       CHARACTER(default_string_length), INTENT(OUT)      :: msg
     582             : 
     583             :       CHARACTER(LEN=20)                                  :: isunit
     584             :       INTEGER                                            :: i, imode, j
     585             : 
     586           1 :       msg = ""
     587           1 :       SELECT CASE (TRIM(myunit))
     588             :       CASE ("eV")
     589           0 :          isunit = "eV"
     590             :       CASE ("K")
     591           1 :          isunit = "K_e"
     592             :       CASE ("atomic energy units ")
     593           0 :          RETURN
     594             :       CASE DEFAULT
     595           0 :          msg = "Unknown unit of C-Matrices for PIGLET. Assuming a.u."
     596           0 :          CPWARN(msg)
     597           1 :          RETURN
     598             :       END SELECT
     599             : 
     600          13 :       DO imode = 1, p
     601         121 :          DO j = 1, nsp1
     602        1092 :             DO i = 1, nsp1
     603        1080 :                c_mat(i, j, imode) = cp_unit_to_cp2k(c_mat(i, j, imode), TRIM(isunit))
     604             :             END DO
     605             :          END DO
     606             :       END DO
     607             : 
     608             :    END SUBROUTINE c_mat_to_cp2k
     609             : 
     610             : ! ***************************************************************************
     611             : !> \brief checks if the matrices are suited for the target temperature
     612             : !> \param line ...
     613             : !> \param propagator ...
     614             : !> \param targettemp ...
     615             : !> \author Felix Uhl
     616             : ! **************************************************************************************************
     617           1 :    SUBROUTINE check_temperature(line, propagator, targettemp)
     618             : 
     619             :       CHARACTER(len=*), INTENT(IN)                       :: line
     620             :       INTEGER, INTENT(IN)                                :: propagator
     621             :       REAL(KIND=dp), INTENT(IN)                          :: targettemp
     622             : 
     623             :       CHARACTER(len=default_string_length)               :: msg
     624             :       INTEGER                                            :: posnumber
     625             :       REAL(KIND=dp)                                      :: convttemp, deviation, matrixtemp, ttemp
     626             : 
     627           1 :       deviation = 100.0d0
     628           1 :       posnumber = INDEX(line, "T=") + 2
     629           1 :       IF (propagator == propagator_rpmd) ttemp = targettemp
     630             :       !Get the matrix temperature
     631           1 :       READ (line(posnumber:), *) matrixtemp
     632           1 :       msg = ""
     633           1 :       IF (INDEX(line, "K") /= 0) THEN
     634           1 :          convttemp = cp_unit_from_cp2k(ttemp, "K")
     635           1 :          IF (ABS(convttemp - matrixtemp) > convttemp/deviation) THEN
     636           0 :             WRITE (UNIT=msg, FMT=*) "PIGLET Simulation temperature (", &
     637           0 :                convttemp, "K) /= matrix temperature (", matrixtemp, "K)"
     638           0 :             CPWARN(msg)
     639             :          END IF
     640           0 :       ELSE IF (INDEX(line, "eV") /= 0) THEN
     641           0 :          convttemp = cp_unit_from_cp2k(ttemp, "K")/11604.505_dp
     642           0 :          IF (ABS(convttemp - matrixtemp) > convttemp/deviation) THEN
     643           0 :             WRITE (UNIT=msg, FMT=*) "PIGLET Simulation temperature (", &
     644           0 :                convttemp, "K) /= matrix temperature (", matrixtemp, "K)"
     645           0 :             CPWARN(msg)
     646             :          END IF
     647           0 :       ELSE IF (INDEX(line, "atomic energy units") /= 0) THEN
     648           0 :          convttemp = ttemp
     649           0 :          IF (ABS(convttemp - matrixtemp) > convttemp/deviation) THEN
     650           0 :             WRITE (UNIT=msg, FMT=*) "PIGLET Simulation temperature (", &
     651           0 :                convttemp, "K) /= matrix temperature (", matrixtemp, "K)"
     652           0 :             CPWARN(msg)
     653             :          END IF
     654             :       ELSE
     655           0 :          WRITE (UNIT=msg, FMT=*) "Unknown PIGLET matrix temperature. Assuming a.u."
     656           0 :          CPWARN(msg)
     657           0 :          convttemp = ttemp
     658           0 :          IF (ABS(convttemp - matrixtemp) > convttemp/deviation) THEN
     659           0 :             WRITE (UNIT=msg, FMT=*) "PIGLET Simulation temperature (", &
     660           0 :                convttemp, "K) /= matrix temperature (", matrixtemp, "K)"
     661           0 :             CPWARN(msg)
     662             :          END IF
     663             :       END IF
     664             : 
     665           1 :    END SUBROUTINE check_temperature
     666             : 
     667             : ! ***************************************************************************
     668             : !> \brief returns the piglet kinetic energy contribution
     669             : !> \param pint_env ...
     670             : !> \author Felix Uhl
     671             : ! **************************************************************************************************
     672          12 :    ELEMENTAL SUBROUTINE pint_calc_piglet_energy(pint_env)
     673             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
     674             : 
     675          12 :       IF (ASSOCIATED(pint_env%piglet_therm)) THEN
     676          12 :          pint_env%e_piglet = pint_env%piglet_therm%thermostat_energy
     677             :       ELSE
     678           0 :          pint_env%e_piglet = 0.0d0
     679             :       END IF
     680             : 
     681          12 :    END SUBROUTINE pint_calc_piglet_energy
     682             : 
     683             : ! ***************************************************************************
     684             : !> \brief calculates S from S*TRANSPOSED(S) by diagonalizing
     685             : !>        if S*TRANSPOSED(S) is a positive definite matrix
     686             : !> \param n order of input matrix
     687             : !> \param SST matrix to be decomposed
     688             : !> \param S result matrix
     689             : !> \author Felix Uhl
     690             : ! **************************************************************************************************
     691          48 :    SUBROUTINE sqrt_pos_def_mat(n, SST, S)
     692             : 
     693             :       INTEGER, INTENT(IN)                                :: n
     694             :       REAL(KIND=dp), DIMENSION(n, n), INTENT(IN)         :: SST
     695             :       REAL(KIND=dp), DIMENSION(n, n), INTENT(OUT)        :: S
     696             : 
     697             :       INTEGER                                            :: i, info, lwork
     698          48 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: work
     699             :       REAL(KIND=dp), DIMENSION(1)                        :: tmplwork
     700          96 :       REAL(KIND=dp), DIMENSION(n)                        :: eigval
     701          48 :       REAL(KIND=dp), DIMENSION(n, n)                     :: A, tmpmatrix
     702             : 
     703             : ! order of input matrix
     704             : ! matrix to be decomposed
     705             : ! result matrix
     706             : ! Variables for eigenvalue/vector computation
     707             : ! store matrix here to pass to lapack routine DSYEVR
     708             : ! Array to contain the eigenvalues
     709             : ! size of temporary real work array
     710             : ! temporary real work array
     711             : ! information about success
     712             : ! counter
     713             : 
     714         480 :       eigval(:) = 0.0_dp
     715        4368 :       A(:, :) = 0.0_dp
     716        4368 :       A(:, :) = SST(:, :)
     717             : 
     718             :       !first call to figure out how big work array needs to be
     719             :       CALL DSYEV('V', & ! Compute eigenvalues and eigenvectors
     720             :                  'U', & ! Store upper triagonal matrix
     721             :                  n, & ! order of matrix A to calculate the eigenvalues/vectors from
     722             :                  A, & ! Matrix to calculate the eigenvalues/vectors from
     723             :                  n, & ! leading order of matrix A
     724             :                  eigval, & ! Array to contain the eigenvalues
     725             :                  tmplwork, & ! temporary real work array
     726             :                  -1, & ! size of temporary real work array
     727          48 :                  info) ! information about success
     728             : 
     729          48 :       lwork = INT(tmplwork(1) + 0.5_dp)
     730         144 :       ALLOCATE (work(lwork))
     731       14736 :       work(:) = 0.0_dp
     732             : 
     733             :       CALL DSYEV('V', & ! Compute eigenvalues and eigenvectors
     734             :                  'U', & ! Store upper triagonal matrix
     735             :                  n, & ! order of matrix A to calculate the eigenvalues/vectors from
     736             :                  A, & ! Matrix to calculate the eigenvalues/vectors from
     737             :                  n, & ! leading order of matrix A
     738             :                  eigval, & ! Array to contain the eigenvalues
     739             :                  work, & ! temporary real work array
     740             :                  lwork, & ! size of temporary real work array
     741          48 :                  info) ! information about success
     742          48 :       DEALLOCATE (work)
     743             :       ! A-matrix now contains the eigenvectors
     744             : 
     745        4368 :       S(:, :) = 0.0_dp
     746         480 :       DO i = 1, n
     747             :          ! In case that numerics made some eigenvalues negative
     748         480 :          IF (eigval(i) > 0.0_dp) THEN
     749         432 :             S(i, i) = SQRT(eigval(i))
     750             :          END IF
     751             :       END DO
     752             : 
     753        4368 :       tmpmatrix(:, :) = 0.0_dp
     754             :       ! Transform matrix back
     755             :       !tmpmatrix = A*S
     756             :       CALL DGEMM('N', & ! A-matrix is not transposed
     757             :                  'N', & ! S-matrix is not transposed
     758             :                  n, & ! number of rows of A-matrix
     759             :                  n, & ! number of columns of S-matrix
     760             :                  n, & ! number of columns of A-matrix and number of rows of S-matrix
     761             :                  1.0D0, & ! scaling factor of A-matrix
     762             :                  A, & ! A-matrix
     763             :                  n, & ! leading dimension of A-matrix as declared
     764             :                  S, & ! S-matrix
     765             :                  n, & ! leading dimension of S-matrix as declared
     766             :                  0.0D0, & ! scaling of tmpmatrix as additive
     767             :                  tmpmatrix, & ! result matrix: tmpmatrix
     768          48 :                  n) ! leading dimension of tmpmatrix
     769             :       !S = tmpmatrix*TRANSPOSE(A) = A*S*TRANSPOSE(A)
     770             :       CALL DGEMM('N', & ! tmpmatrix not transposed
     771             :                  'T', & ! A-matrix is transposed
     772             :                  n, & ! number of rows of tmpmatrix
     773             :                  n, & ! number of columns of A-matrix
     774             :                  n, & ! number of columns of tmpmatrix and rows of A-matrix
     775             :                  1.0D0, & ! scaling factor of tmpmatrix
     776             :                  tmpmatrix, & ! tmpmatrix
     777             :                  n, & ! leading dimension of tmpmatrix as declared
     778             :                  A, & ! A-matrix
     779             :                  n, & ! leading dimension of A-matrix as declared
     780             :                  0.0D0, & ! scaling of S-matrix as additive
     781             :                  S, & ! result matrix: S-matrix
     782          48 :                  n) ! leading dimension of S-matrix as declared
     783             : 
     784          48 :    END SUBROUTINE sqrt_pos_def_mat
     785             : 
     786             : END MODULE pint_piglet

Generated by: LCOV version 1.15