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 Methods to apply GLE to PI runs.
10 : !> \author michelec
11 : !> \par History
12 : !> 06.2010 created [michelec]
13 : !> \note trying to keep duplication at a minimum....
14 : ! **************************************************************************************************
15 :
16 : MODULE pint_gle
17 : USE gle_system_dynamics, ONLY: gle_cholesky_stab
18 : USE gle_system_types, ONLY: gle_type
19 : USE kinds, ONLY: dp
20 : USE pint_types, ONLY: pint_env_type
21 : #include "../base/base_uses.f90"
22 :
23 : IMPLICIT NONE
24 :
25 : PRIVATE
26 :
27 : PUBLIC :: pint_gle_step, pint_gle_init, pint_calc_gle_energy
28 :
29 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pint_gle'
30 :
31 : CONTAINS
32 : ! **************************************************************************************************
33 : !> \brief ...
34 : !> \param pint_env ...
35 : ! **************************************************************************************************
36 6 : ELEMENTAL SUBROUTINE pint_calc_gle_energy(pint_env)
37 : TYPE(pint_env_type), INTENT(INOUT) :: pint_env
38 :
39 : INTEGER :: i
40 :
41 6 : pint_env%e_gle = 0._dp
42 6 : IF (ASSOCIATED(pint_env%gle)) THEN
43 55302 : DO i = 1, pint_env%gle%loc_num_gle
44 55302 : pint_env%e_gle = pint_env%e_gle + pint_env%gle%nvt(i)%thermostat_energy
45 : END DO
46 : END IF
47 6 : END SUBROUTINE
48 :
49 : ! **************************************************************************************************
50 : !> \brief ...
51 : !> \param pint_env ...
52 : ! **************************************************************************************************
53 2 : SUBROUTINE pint_gle_init(pint_env)
54 : TYPE(pint_env_type), INTENT(INOUT) :: pint_env
55 :
56 : INTEGER :: i, ib, idim, imap, j
57 4 : REAL(dp) :: mf, rr(pint_env%gle%ndim), cc(pint_env%gle%ndim, pint_env%gle%ndim)
58 :
59 2 : CALL gle_cholesky_stab(pint_env%gle%c_mat, cc, pint_env%gle%ndim)
60 18434 : DO i = 1, pint_env%gle%loc_num_gle
61 18432 : imap = pint_env%gle%map_info%index(i)
62 18432 : ib = 1 + (imap - 1)/pint_env%ndim
63 18432 : idim = 1 + MOD(imap - 1, pint_env%ndim)
64 18432 : mf = 1.0_dp/SQRT(pint_env%mass_fict(ib, idim))
65 110592 : DO j = 1, pint_env%gle%ndim
66 110592 : rr(j) = pint_env%gle%nvt(i)%gaussian_rng_stream%next()*mf
67 : END DO
68 663554 : pint_env%gle%nvt(i)%s = MATMUL(cc, rr)
69 : END DO
70 :
71 2 : END SUBROUTINE pint_gle_init
72 :
73 : ! **************************************************************************************************
74 : !> \brief ...
75 : !> \param pint_env ...
76 : ! **************************************************************************************************
77 4 : SUBROUTINE pint_gle_step(pint_env)
78 : TYPE(pint_env_type), INTENT(INOUT) :: pint_env
79 :
80 : CHARACTER(len=*), PARAMETER :: routineN = 'pint_gle_step'
81 :
82 : INTEGER :: handle, iadd, ib, ideg, idim, imap, &
83 : ndim, num
84 : REAL(dp) :: alpha, beta, mf, rr
85 4 : REAL(dp), DIMENSION(:, :), POINTER :: a_mat, e_tmp, h_tmp, s_tmp
86 : TYPE(gle_type), POINTER :: gle
87 :
88 4 : CALL timeset(routineN, handle)
89 :
90 4 : gle => pint_env%gle
91 4 : ndim = gle%ndim
92 :
93 16 : ALLOCATE (s_tmp(ndim, gle%loc_num_gle))
94 221188 : s_tmp = 0.0_dp
95 16 : ALLOCATE (e_tmp(ndim, gle%loc_num_gle))
96 16 : ALLOCATE (h_tmp(ndim, gle%loc_num_gle))
97 :
98 36868 : DO ideg = 1, gle%loc_num_gle
99 36864 : imap = gle%map_info%index(ideg)
100 36864 : ib = 1 + (imap - 1)/pint_env%ndim
101 36864 : idim = 1 + MOD(imap - 1, pint_env%ndim)
102 :
103 36864 : gle%nvt(ideg)%s(1) = pint_env%uv_t(ib, idim)
104 : gle%nvt(ideg)%thermostat_energy = gle%nvt(ideg)%thermostat_energy &
105 36864 : + 0.5_dp*pint_env%mass_fict(ib, idim)*gle%nvt(ideg)%s(1)**2
106 36864 : s_tmp(1, imap) = gle%nvt(ideg)%s(1)
107 36864 : rr = gle%nvt(ideg)%gaussian_rng_stream%next()
108 36864 : mf = 1.0_dp/SQRT(pint_env%mass_fict(ib, idim))
109 36864 : e_tmp(1, imap) = rr*mf
110 184324 : DO iadd = 2, ndim
111 147456 : s_tmp(iadd, imap) = gle%nvt(ideg)%s(iadd)
112 147456 : rr = gle%nvt(ideg)%gaussian_rng_stream%next()
113 184320 : e_tmp(iadd, imap) = rr*mf
114 : END DO
115 : END DO
116 4 : num = gle%loc_num_gle
117 4 : a_mat => gle%gle_s
118 4 : alpha = 1.0_dp
119 4 : beta = 0.0_dp
120 :
121 4 : CALL DGEMM('N', 'N', ndim, num, ndim, alpha, a_mat(1, 1), ndim, e_tmp(1, 1), ndim, beta, h_tmp(1, 1), ndim)
122 :
123 4 : a_mat => gle%gle_t
124 4 : beta = 1.0_dp
125 4 : CALL DGEMM("N", "N", ndim, num, ndim, alpha, a_mat(1, 1), ndim, s_tmp(1, 1), ndim, beta, h_tmp(1, 1), ndim)
126 :
127 36868 : DO ideg = 1, gle%loc_num_gle
128 36864 : imap = gle%map_info%index(ideg)
129 :
130 221184 : DO iadd = 1, ndim
131 221184 : gle%nvt(ideg)%s(iadd) = h_tmp(iadd, imap)
132 : END DO
133 36864 : ib = 1 + (imap - 1)/pint_env%ndim
134 36864 : idim = 1 + MOD(imap - 1, pint_env%ndim)
135 36864 : pint_env%uv_t(ib, idim) = gle%nvt(ideg)%s(1)
136 : gle%nvt(ideg)%thermostat_energy = gle%nvt(ideg)%thermostat_energy &
137 36868 : - 0.5_dp*pint_env%mass_fict(ib, idim)*gle%nvt(ideg)%s(1)**2
138 : END DO
139 4 : pint_env%e_kin_t = 0.0_dp
140 4 : DEALLOCATE (e_tmp, s_tmp, h_tmp)
141 4 : CALL timestop(handle)
142 4 : END SUBROUTINE
143 : END MODULE
|
