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 A collection of methods to treat the QM/MM links
10 : !> \par History
11 : !> 12.2004 created [tlaino]
12 : !> \author Teodoro Laino
13 : ! **************************************************************************************************
14 : MODULE qmmm_links_methods
15 :
16 : USE cp_log_handling, ONLY: cp_to_string
17 : USE kinds, ONLY: dp
18 : USE particle_types, ONLY: particle_type
19 : USE qmmm_types_low, ONLY: add_set_type,&
20 : qmmm_env_qm_type,&
21 : qmmm_imomm_link_type,&
22 : qmmm_links_type
23 : #include "./base/base_uses.f90"
24 :
25 : IMPLICIT NONE
26 : PRIVATE
27 :
28 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
29 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_links_methods'
30 : PUBLIC :: qmmm_link_Imomm_coord, &
31 : qmmm_link_Imomm_forces, &
32 : qmmm_added_chrg_coord, &
33 : qmmm_added_chrg_forces
34 :
35 : CONTAINS
36 :
37 : ! **************************************************************************************************
38 : !> \brief correct the position for qm/mm IMOMM link type
39 : !> \param qmmm_links ...
40 : !> \param particles ...
41 : !> \param qm_atom_index ...
42 : !> \par History
43 : !> 12.2004 created [tlaino]
44 : !> \author Teodoro Laino
45 : ! **************************************************************************************************
46 338 : SUBROUTINE qmmm_link_Imomm_coord(qmmm_links, particles, qm_atom_index)
47 : TYPE(qmmm_links_type), POINTER :: qmmm_links
48 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
49 : INTEGER, DIMENSION(:), POINTER :: qm_atom_index
50 :
51 : INTEGER :: ilink, ip, ip_mm, ip_qm, mm_index, &
52 : n_imomm, qm_index
53 : REAL(KIND=dp) :: alpha
54 : TYPE(qmmm_imomm_link_type), POINTER :: my_link
55 :
56 338 : n_imomm = SIZE(qmmm_links%imomm)
57 338 : CPASSERT(n_imomm /= 0)
58 1270 : DO ilink = 1, n_imomm
59 932 : my_link => qmmm_links%imomm(ilink)%link
60 932 : qm_index = my_link%qm_index
61 932 : mm_index = my_link%mm_index
62 932 : alpha = 1.0_dp/my_link%alpha
63 25616 : DO ip = 1, SIZE(qm_atom_index)
64 25616 : IF (qm_atom_index(ip) == qm_index) EXIT
65 : END DO
66 932 : IF (ip == SIZE(qm_atom_index) + 1) &
67 : CALL cp_abort(__LOCATION__, &
68 : "QM atom index ("//cp_to_string(qm_index)//") specified in the LINK section nr.("// &
69 0 : cp_to_string(ilink)//") is not defined as a QM atom! Please inspect your QM_KIND sections. ")
70 932 : ip_qm = ip
71 44324 : DO ip = 1, SIZE(qm_atom_index)
72 44324 : IF (qm_atom_index(ip) == mm_index) EXIT
73 : END DO
74 932 : IF (ip == SIZE(qm_atom_index) + 1) &
75 : CALL cp_abort(__LOCATION__, &
76 : "Error in setting up the MM atom index ("//cp_to_string(mm_index)// &
77 0 : ") specified in the LINK section nr.("//cp_to_string(ilink)//"). Please report this bug! ")
78 932 : ip_mm = ip
79 4066 : particles(ip_mm)%r = alpha*particles(ip_mm)%r + (1.0_dp - alpha)*particles(ip_qm)%r
80 : END DO
81 :
82 338 : END SUBROUTINE qmmm_link_Imomm_coord
83 :
84 : ! **************************************************************************************************
85 : !> \brief correct the forces for qm/mm IMOMM link type
86 : !> \param qmmm_links ...
87 : !> \param particles_qm ...
88 : !> \param qm_atom_index ...
89 : !> \par History
90 : !> 12.2004 created [tlaino]
91 : !> \author Teodoro Laino
92 : ! **************************************************************************************************
93 274 : SUBROUTINE qmmm_link_Imomm_forces(qmmm_links, particles_qm, qm_atom_index)
94 : TYPE(qmmm_links_type), POINTER :: qmmm_links
95 : TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm
96 : INTEGER, DIMENSION(:), POINTER :: qm_atom_index
97 :
98 : INTEGER :: ilink, ip, ip_mm, ip_qm, mm_index, &
99 : n_imomm, qm_index
100 : REAL(KIND=dp) :: alpha
101 : TYPE(qmmm_imomm_link_type), POINTER :: my_link
102 :
103 274 : n_imomm = SIZE(qmmm_links%imomm)
104 274 : CPASSERT(n_imomm /= 0)
105 1008 : DO ilink = 1, n_imomm
106 734 : my_link => qmmm_links%imomm(ilink)%link
107 734 : qm_index = my_link%qm_index
108 734 : mm_index = my_link%mm_index
109 734 : alpha = 1.0_dp/my_link%alpha
110 19558 : DO ip = 1, SIZE(qm_atom_index)
111 19558 : IF (qm_atom_index(ip) == qm_index) EXIT
112 : END DO
113 734 : IF (ip == SIZE(qm_atom_index) + 1) &
114 : CALL cp_abort(__LOCATION__, &
115 : "QM atom index ("//cp_to_string(qm_index)//") specified in the LINK section nr.("// &
116 0 : cp_to_string(ilink)//") is not defined as a QM atom! Please inspect your QM_KIND sections. ")
117 734 : ip_qm = ip
118 34278 : DO ip = 1, SIZE(qm_atom_index)
119 34278 : IF (qm_atom_index(ip) == mm_index) EXIT
120 : END DO
121 734 : IF (ip == SIZE(qm_atom_index) + 1) &
122 : CALL cp_abort(__LOCATION__, &
123 : "Error in setting up the MM atom index ("//cp_to_string(mm_index)// &
124 0 : ") specified in the LINK section nr.("//cp_to_string(ilink)//"). Please report this bug! ")
125 734 : ip_mm = ip
126 2936 : particles_qm(ip_qm)%f = particles_qm(ip_qm)%f + particles_qm(ip_mm)%f*(1.0_dp - alpha)
127 3210 : particles_qm(ip_mm)%f = particles_qm(ip_mm)%f*alpha
128 : END DO
129 :
130 274 : END SUBROUTINE qmmm_link_Imomm_forces
131 :
132 : ! **************************************************************************************************
133 : !> \brief correct the position for added charges in qm/mm link scheme
134 : !> \param qmmm_env ...
135 : !> \param particles ...
136 : !> \par History
137 : !> 01.2005 created [tlaino]
138 : !> \author Teodoro Laino
139 : ! **************************************************************************************************
140 32 : SUBROUTINE qmmm_added_chrg_coord(qmmm_env, particles)
141 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
142 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
143 :
144 : INTEGER :: I, Index1, Index2
145 : REAL(KIND=dp) :: alpha
146 : TYPE(add_set_type), POINTER :: added_charges
147 :
148 32 : added_charges => qmmm_env%added_charges
149 :
150 144 : DO i = 1, added_charges%num_mm_atoms
151 112 : Index1 = added_charges%add_env(i)%Index1
152 112 : Index2 = added_charges%add_env(i)%Index2
153 112 : alpha = added_charges%add_env(i)%alpha
154 928 : added_charges%added_particles(i)%r = alpha*particles(Index1)%r + (1.0_dp - alpha)*particles(Index2)%r
155 : END DO
156 :
157 32 : END SUBROUTINE qmmm_added_chrg_coord
158 :
159 : ! **************************************************************************************************
160 : !> \brief correct the forces due to the added charges in qm/mm link scheme
161 : !> \param qmmm_env ...
162 : !> \param particles ...
163 : !> \par History
164 : !> 01.2005 created [tlaino]
165 : !> \author Teodoro Laino
166 : ! **************************************************************************************************
167 32 : SUBROUTINE qmmm_added_chrg_forces(qmmm_env, particles)
168 : TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
169 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
170 :
171 : INTEGER :: I, Index1, Index2
172 : REAL(KIND=dp) :: alpha
173 : TYPE(add_set_type), POINTER :: added_charges
174 :
175 32 : added_charges => qmmm_env%added_charges
176 :
177 144 : DO i = 1, added_charges%num_mm_atoms
178 112 : Index1 = added_charges%add_env(i)%Index1
179 112 : Index2 = added_charges%add_env(i)%Index2
180 112 : alpha = added_charges%add_env(i)%alpha
181 896 : particles(Index1)%f = particles(Index1)%f + alpha*added_charges%added_particles(i)%f
182 928 : particles(Index2)%f = particles(Index2)%f + (1.0_dp - alpha)*added_charges%added_particles(i)%f
183 : END DO
184 :
185 32 : END SUBROUTINE qmmm_added_chrg_forces
186 :
187 : END MODULE qmmm_links_methods
|
