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 used in the flexible partitioning scheme
10 : !> \par History
11 : !> 04.2006 [Joost VandeVondele]
12 : !> \author Joost VandeVondele
13 : ! **************************************************************************************************
14 : MODULE fp_methods
15 :
16 : USE beta_gamma_psi, ONLY: psi
17 : USE cell_types, ONLY: cell_type,&
18 : pbc
19 : USE cp_log_handling, ONLY: cp_get_default_logger,&
20 : cp_logger_type
21 : USE cp_output_handling, ONLY: cp_iter_string,&
22 : cp_print_key_finished_output,&
23 : cp_print_key_unit_nr
24 : USE cp_subsys_types, ONLY: cp_subsys_get,&
25 : cp_subsys_type
26 : USE fp_types, ONLY: fp_type
27 : USE kinds, ONLY: dp
28 : USE mathconstants, ONLY: fac,&
29 : maxfac,&
30 : oorootpi
31 : USE particle_list_types, ONLY: particle_list_type
32 : USE particle_types, ONLY: particle_type
33 : #include "./base/base_uses.f90"
34 :
35 : IMPLICIT NONE
36 : PRIVATE
37 :
38 : PUBLIC :: fp_eval
39 :
40 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fp_methods'
41 :
42 : CONTAINS
43 :
44 : ! **************************************************************************************************
45 : !> \brief computest the forces and the energy due to the flexible potential & bias,
46 : !> and writes the weights file
47 : !> \param fp_env ...
48 : !> \param subsys ...
49 : !> \param cell ...
50 : !> \par History
51 : !> 04.2006 created [Joost VandeVondele]
52 : ! **************************************************************************************************
53 122 : SUBROUTINE fp_eval(fp_env, subsys, cell)
54 : TYPE(fp_type), POINTER :: fp_env
55 : TYPE(cp_subsys_type), POINTER :: subsys
56 : TYPE(cell_type), POINTER :: cell
57 :
58 : CHARACTER(len=*), PARAMETER :: routineN = 'fp_eval'
59 :
60 : CHARACTER(LEN=15) :: tmpstr
61 : INTEGER :: handle, i, icenter, iparticle, &
62 : output_unit
63 : LOGICAL :: zero_weight
64 : REAL(KIND=dp) :: c, dcdr, kT, r, rab(3), sf, strength
65 : TYPE(cp_logger_type), POINTER :: logger
66 : TYPE(particle_list_type), POINTER :: particles_list
67 122 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
68 :
69 122 : CALL timeset(routineN, handle)
70 :
71 122 : CPASSERT(ASSOCIATED(fp_env))
72 122 : CPASSERT(fp_env%use_fp)
73 122 : CPASSERT(ASSOCIATED(subsys))
74 122 : CALL cp_subsys_get(subsys, particles=particles_list)
75 122 : particles => particles_list%els
76 :
77 : ! compute the force due to the reflecting walls
78 : ! and count the distribution in discrete and contiguous ways
79 122 : zero_weight = .FALSE.
80 122 : fp_env%restraint_energy = 0.0_dp
81 122 : icenter = fp_env%central_atom
82 122 : strength = fp_env%strength
83 122 : fp_env%i1 = 0; fp_env%i2 = 0; fp_env%o1 = 0; fp_env%o2 = 0
84 122 : fp_env%ri1 = 0.0_dp; fp_env%ri2 = 0.0_dp; fp_env%ro1 = 0.0_dp; fp_env%ro2 = 0.0_dp
85 122 : fp_env%energy = 0.0_dp
86 :
87 : ! inner particles
88 1098 : DO i = 1, SIZE(fp_env%inner_atoms)
89 976 : iparticle = fp_env%inner_atoms(i)
90 3904 : rab = particles(iparticle)%r - particles(icenter)%r
91 3904 : rab = pbc(rab, cell)
92 3904 : r = SQRT(SUM(rab**2))
93 : ! constraint wall (they feel to outer wall)
94 976 : IF (r > fp_env%outer_radius) THEN
95 0 : zero_weight = .TRUE.
96 0 : fp_env%restraint_energy = fp_env%restraint_energy + 0.5_dp*strength*(r - fp_env%outer_radius)**2
97 0 : sf = strength*(r - fp_env%outer_radius)/r
98 0 : particles(iparticle)%f = particles(iparticle)%f - sf*rab
99 0 : particles(icenter)%f = particles(icenter)%f + sf*rab
100 : END IF
101 : ! count the distribution
102 976 : IF (r > fp_env%inner_radius) THEN
103 610 : fp_env%i2 = fp_env%i2 + 1
104 : ELSE
105 366 : fp_env%i1 = fp_env%i1 + 1
106 : END IF
107 : ! smooth count the distribution
108 976 : CALL smooth_count(r, fp_env%inner_radius, fp_env%smooth_width, c, dcdr)
109 976 : fp_env%ri1 = fp_env%ri1 + c
110 1098 : fp_env%ri2 = fp_env%ri2 + (1.0_dp - c)
111 : END DO
112 :
113 : ! outer particles
114 2928 : DO i = 1, SIZE(fp_env%outer_atoms)
115 2806 : iparticle = fp_env%outer_atoms(i)
116 11224 : rab = particles(iparticle)%r - particles(icenter)%r
117 11224 : rab = pbc(rab, cell)
118 11224 : r = SQRT(SUM(rab**2))
119 : ! constraint wall (they feel the inner wall)
120 2806 : IF (r < fp_env%inner_radius) THEN
121 0 : zero_weight = .TRUE.
122 : fp_env%restraint_energy = fp_env%restraint_energy + &
123 0 : 0.5_dp*strength*(r - fp_env%inner_radius)**2
124 0 : sf = strength*(r - fp_env%inner_radius)/r
125 0 : particles(iparticle)%f = particles(iparticle)%f - sf*rab
126 0 : particles(icenter)%f = particles(icenter)%f + sf*rab
127 : END IF
128 : ! count the distribution
129 2806 : IF (r > fp_env%outer_radius) THEN
130 1998 : fp_env%o2 = fp_env%o2 + 1
131 : ELSE
132 808 : fp_env%o1 = fp_env%o1 + 1
133 : END IF
134 : ! smooth count the distribution
135 2806 : CALL smooth_count(r, fp_env%outer_radius, fp_env%smooth_width, c, dcdr)
136 2806 : fp_env%ro1 = fp_env%ro1 + c
137 2928 : fp_env%ro2 = fp_env%ro2 + (1.0_dp - c)
138 : END DO
139 122 : fp_env%energy = fp_env%energy + fp_env%restraint_energy
140 :
141 : ! the combinatorial weight
142 122 : i = fp_env%i2 + fp_env%o1
143 122 : CPASSERT(i <= maxfac)
144 122 : fp_env%comb_weight = (fac(fp_env%i2)*fac(fp_env%o1))/fac(i)
145 :
146 : ! we can add the bias potential now.
147 : ! this bias has the form
148 : ! kT * { ln[(o1+i2)!] - ln[o1!] - ln[i2!] }
149 : ! where the smooth counts are used for o1 and i2
150 122 : fp_env%bias_energy = 0.0_dp
151 122 : IF (fp_env%bias) THEN
152 122 : kT = fp_env%temperature
153 : fp_env%bias_energy = kT*(LOG_GAMMA(fp_env%ro1 + fp_env%ri2 + 1) - &
154 122 : LOG_GAMMA(fp_env%ro1 + 1) - LOG_GAMMA(fp_env%ri2 + 1))
155 :
156 : ! and add the corresponding forces
157 : ! inner particles
158 1098 : DO i = 1, SIZE(fp_env%inner_atoms)
159 976 : iparticle = fp_env%inner_atoms(i)
160 3904 : rab = particles(iparticle)%r - particles(icenter)%r
161 3904 : rab = pbc(rab, cell)
162 3904 : r = SQRT(SUM(rab**2))
163 976 : CALL smooth_count(r, fp_env%inner_radius, fp_env%smooth_width, c, dcdr)
164 976 : sf = kT*(psi(fp_env%ro1 + fp_env%ri2 + 1) - psi(fp_env%ri2 + 1))*(-dcdr)/r
165 3904 : particles(iparticle)%f = particles(iparticle)%f - sf*rab
166 5002 : particles(icenter)%f = particles(icenter)%f + sf*rab
167 : END DO
168 : ! outer particles
169 2928 : DO i = 1, SIZE(fp_env%outer_atoms)
170 2806 : iparticle = fp_env%outer_atoms(i)
171 11224 : rab = particles(iparticle)%r - particles(icenter)%r
172 11224 : rab = pbc(rab, cell)
173 11224 : r = SQRT(SUM(rab**2))
174 2806 : CALL smooth_count(r, fp_env%outer_radius, fp_env%smooth_width, c, dcdr)
175 2806 : sf = kT*(psi(fp_env%ro1 + fp_env%ri2 + 1) - psi(fp_env%ro1 + 1))*(dcdr)/r
176 11224 : particles(iparticle)%f = particles(iparticle)%f - sf*rab
177 14152 : particles(icenter)%f = particles(icenter)%f + sf*rab
178 : END DO
179 : END IF
180 122 : fp_env%energy = fp_env%energy + fp_env%bias_energy
181 122 : fp_env%bias_weight = EXP(fp_env%bias_energy/kT)
182 :
183 : ! if this configuration is a valid one, compute its weight
184 122 : IF (zero_weight) THEN
185 0 : fp_env%weight = 0.0_dp
186 : ELSE
187 122 : fp_env%weight = fp_env%comb_weight*fp_env%bias_weight
188 : END IF
189 :
190 : ! put weights and other info on file
191 122 : logger => cp_get_default_logger()
192 : output_unit = cp_print_key_unit_nr(logger, fp_env%print_section, "", &
193 122 : extension=".weights")
194 122 : IF (output_unit > 0) THEN
195 61 : tmpstr = cp_iter_string(logger%iter_info, fp_env%print_section)
196 : WRITE (output_unit, '(T2,A15,6(1X,F16.10),4(1X,I4),4(1X,F16.10))') &
197 61 : tmpstr, &
198 61 : fp_env%weight, fp_env%comb_weight, fp_env%bias_weight, &
199 61 : fp_env%energy, fp_env%restraint_energy, fp_env%bias_energy, &
200 61 : fp_env%i1, fp_env%i2, fp_env%o1, fp_env%o2, &
201 122 : fp_env%ri1, fp_env%ri2, fp_env%ro1, fp_env%ro2
202 : END IF
203 :
204 : CALL cp_print_key_finished_output(output_unit, logger, fp_env%print_section, &
205 122 : "")
206 :
207 122 : CALL timestop(handle)
208 :
209 122 : END SUBROUTINE fp_eval
210 :
211 : ! **************************************************************************************************
212 : !> \brief counts in a smooth way (error function with width=width)
213 : !> if r is closer than r1. Returns 1.0 for the count=c if r<<r1
214 : !> and the derivative wrt r dcdr
215 : !> \param r ...
216 : !> \param r1 ...
217 : !> \param width ...
218 : !> \param c ...
219 : !> \param dcdr ...
220 : !> \par History
221 : !> 04.2006 created [Joost VandeVondele]
222 : ! **************************************************************************************************
223 7564 : SUBROUTINE smooth_count(r, r1, width, c, dcdr)
224 : REAL(KIND=dp), INTENT(IN) :: r, r1, width
225 : REAL(KIND=dp), INTENT(OUT) :: c, dcdr
226 :
227 : REAL(KIND=dp) :: arg
228 :
229 7564 : arg = (r1 - r)/width
230 :
231 7564 : c = (1.0_dp + ERF(arg))/2.0_dp
232 7564 : dcdr = (-oorootpi/width)*EXP(-arg**2)
233 :
234 7564 : END SUBROUTINE
235 :
236 : END MODULE fp_methods
|
