(git:e7e05ae)
fp_methods.F
Go to the documentation of this file.
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  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  TYPE(particle_type), DIMENSION(:), POINTER :: particles
68 
69  CALL timeset(routinen, handle)
70 
71  cpassert(ASSOCIATED(fp_env))
72  cpassert(fp_env%use_fp)
73  cpassert(ASSOCIATED(subsys))
74  CALL cp_subsys_get(subsys, particles=particles_list)
75  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  zero_weight = .false.
80  fp_env%restraint_energy = 0.0_dp
81  icenter = fp_env%central_atom
82  strength = fp_env%strength
83  fp_env%i1 = 0; fp_env%i2 = 0; fp_env%o1 = 0; fp_env%o2 = 0
84  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  fp_env%energy = 0.0_dp
86 
87  ! inner particles
88  DO i = 1, SIZE(fp_env%inner_atoms)
89  iparticle = fp_env%inner_atoms(i)
90  rab = particles(iparticle)%r - particles(icenter)%r
91  rab = pbc(rab, cell)
92  r = sqrt(sum(rab**2))
93  ! constraint wall (they feel to outer wall)
94  IF (r > fp_env%outer_radius) THEN
95  zero_weight = .true.
96  fp_env%restraint_energy = fp_env%restraint_energy + 0.5_dp*strength*(r - fp_env%outer_radius)**2
97  sf = strength*(r - fp_env%outer_radius)/r
98  particles(iparticle)%f = particles(iparticle)%f - sf*rab
99  particles(icenter)%f = particles(icenter)%f + sf*rab
100  END IF
101  ! count the distribution
102  IF (r > fp_env%inner_radius) THEN
103  fp_env%i2 = fp_env%i2 + 1
104  ELSE
105  fp_env%i1 = fp_env%i1 + 1
106  END IF
107  ! smooth count the distribution
108  CALL smooth_count(r, fp_env%inner_radius, fp_env%smooth_width, c, dcdr)
109  fp_env%ri1 = fp_env%ri1 + c
110  fp_env%ri2 = fp_env%ri2 + (1.0_dp - c)
111  END DO
112 
113  ! outer particles
114  DO i = 1, SIZE(fp_env%outer_atoms)
115  iparticle = fp_env%outer_atoms(i)
116  rab = particles(iparticle)%r - particles(icenter)%r
117  rab = pbc(rab, cell)
118  r = sqrt(sum(rab**2))
119  ! constraint wall (they feel the inner wall)
120  IF (r < fp_env%inner_radius) THEN
121  zero_weight = .true.
122  fp_env%restraint_energy = fp_env%restraint_energy + &
123  0.5_dp*strength*(r - fp_env%inner_radius)**2
124  sf = strength*(r - fp_env%inner_radius)/r
125  particles(iparticle)%f = particles(iparticle)%f - sf*rab
126  particles(icenter)%f = particles(icenter)%f + sf*rab
127  END IF
128  ! count the distribution
129  IF (r > fp_env%outer_radius) THEN
130  fp_env%o2 = fp_env%o2 + 1
131  ELSE
132  fp_env%o1 = fp_env%o1 + 1
133  END IF
134  ! smooth count the distribution
135  CALL smooth_count(r, fp_env%outer_radius, fp_env%smooth_width, c, dcdr)
136  fp_env%ro1 = fp_env%ro1 + c
137  fp_env%ro2 = fp_env%ro2 + (1.0_dp - c)
138  END DO
139  fp_env%energy = fp_env%energy + fp_env%restraint_energy
140 
141  ! the combinatorial weight
142  i = fp_env%i2 + fp_env%o1
143  cpassert(i <= maxfac)
144  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  fp_env%bias_energy = 0.0_dp
151  IF (fp_env%bias) THEN
152  kt = fp_env%temperature
153  fp_env%bias_energy = kt*(log_gamma(fp_env%ro1 + fp_env%ri2 + 1) - &
154  log_gamma(fp_env%ro1 + 1) - log_gamma(fp_env%ri2 + 1))
155 
156  ! and add the corresponding forces
157  ! inner particles
158  DO i = 1, SIZE(fp_env%inner_atoms)
159  iparticle = fp_env%inner_atoms(i)
160  rab = particles(iparticle)%r - particles(icenter)%r
161  rab = pbc(rab, cell)
162  r = sqrt(sum(rab**2))
163  CALL smooth_count(r, fp_env%inner_radius, fp_env%smooth_width, c, dcdr)
164  sf = kt*(psi(fp_env%ro1 + fp_env%ri2 + 1) - psi(fp_env%ri2 + 1))*(-dcdr)/r
165  particles(iparticle)%f = particles(iparticle)%f - sf*rab
166  particles(icenter)%f = particles(icenter)%f + sf*rab
167  END DO
168  ! outer particles
169  DO i = 1, SIZE(fp_env%outer_atoms)
170  iparticle = fp_env%outer_atoms(i)
171  rab = particles(iparticle)%r - particles(icenter)%r
172  rab = pbc(rab, cell)
173  r = sqrt(sum(rab**2))
174  CALL smooth_count(r, fp_env%outer_radius, fp_env%smooth_width, c, dcdr)
175  sf = kt*(psi(fp_env%ro1 + fp_env%ri2 + 1) - psi(fp_env%ro1 + 1))*(dcdr)/r
176  particles(iparticle)%f = particles(iparticle)%f - sf*rab
177  particles(icenter)%f = particles(icenter)%f + sf*rab
178  END DO
179  END IF
180  fp_env%energy = fp_env%energy + fp_env%bias_energy
181  fp_env%bias_weight = exp(fp_env%bias_energy/kt)
182 
183  ! if this configuration is a valid one, compute its weight
184  IF (zero_weight) THEN
185  fp_env%weight = 0.0_dp
186  ELSE
187  fp_env%weight = fp_env%comb_weight*fp_env%bias_weight
188  END IF
189 
190  ! put weights and other info on file
191  logger => cp_get_default_logger()
192  output_unit = cp_print_key_unit_nr(logger, fp_env%print_section, "", &
193  extension=".weights")
194  IF (output_unit > 0) THEN
195  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  tmpstr, &
198  fp_env%weight, fp_env%comb_weight, fp_env%bias_weight, &
199  fp_env%energy, fp_env%restraint_energy, fp_env%bias_energy, &
200  fp_env%i1, fp_env%i2, fp_env%o1, fp_env%o2, &
201  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  "")
206 
207  CALL timestop(handle)
208 
209  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  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  arg = (r1 - r)/width
230 
231  c = (1.0_dp + erf(arg))/2.0_dp
232  dcdr = (-oorootpi/width)*exp(-arg**2)
233 
234  END SUBROUTINE
235 
236 END MODULE fp_methods
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
beta_gamma_psi
Definition: beta_gamma_psi.F:9
beta_gamma_psi::psi
real(dp) function, public psi(xx)
...
Definition: beta_gamma_psi.F:718
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_iter_string
character(len=default_string_length) function, public cp_iter_string(iter_info, print_key, for_file)
returns the iteration string, a string that is useful to create unique filenames (once you trim it)
Definition: cp_output_handling.F:517
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition: cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
Definition: cp_subsys_types.F:345
fp_methods
methods used in the flexible partitioning scheme
Definition: fp_methods.F:14
fp_methods::fp_eval
subroutine, public fp_eval(fp_env, subsys, cell)
computest the forces and the energy due to the flexible potential & bias, and writes the weights file
Definition: fp_methods.F:54
fp_types
types used in the flexible partitioning scheme
Definition: fp_types.F:14
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::oorootpi
real(kind=dp), parameter, public oorootpi
Definition: mathconstants.F:115
mathconstants::maxfac
integer, parameter, public maxfac
Definition: mathconstants.F:36
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition: mathconstants.F:37
particle_list_types
represent a simple array based list of the given type
Definition: particle_list_types.F:15
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19