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