1! ------------ ---------------------------------------------------------- 2! LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator 3! https://www.lammps.org/ Sandia National Laboratories 4! Steve Plimpton, sjplimp@sandia.gov 5! 6! Copyright (2003) Sandia Corporation. Under the terms of Contract 7! DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains 8! certain rights in this software. This software is distributed under 9! the GNU General Public License. 10! 11! See the README file in the top-level LAMMPS directory. 12! 13! Contributing author: Alexey N. Volkov, UA, avolkov1@ua.edu 14!------------------------------------------------------------------------- 15 16module TPMM1 !************************************************************************************** 17! 18! Combined/Weighted potential of type 1. 19! 20! Calculation of the combined potential is based on the 'extended' chain. 21! 22!--------------------------------------------------------------------------------------------------- 23! 24! Intel Fortran. 25! 26! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 09.01, 2017 27! 28!*************************************************************************************************** 29 30use TubePotMono 31use iso_c_binding, only : c_int, c_double, c_char 32implicit none 33 34!--------------------------------------------------------------------------------------------------- 35! Constants 36!--------------------------------------------------------------------------------------------------- 37 38 ! Maximum length of a segment chain 39 integer(c_int), parameter :: TPM_MAX_CHAIN = 100 40 41!--------------------------------------------------------------------------------------------------- 42! Numerical parameters 43!--------------------------------------------------------------------------------------------------- 44 45 ! Switching parameters 46 real(c_double) :: TPMC123 = 1.0d+00 ! Non-dimensional 47 real(c_double) :: TPMC3 = 10.0d+00 ! (A) 48 49!--------------------------------------------------------------------------------------------------- 50! Global variables 51!--------------------------------------------------------------------------------------------------- 52 53 ! These global variables are used to speedup calculations 54 real(c_double), dimension(0:2,0:TPM_MAX_CHAIN-1) :: E1, E2, EE1, EE2 55 real(c_double), dimension(0:2) :: Q1, Q2, Qe, Qe1, DR, Z1, Z2, S1, S2, Pe, Pe1 56 real(c_double), dimension(0:TPM_MAX_CHAIN-1) :: W, C 57 real(c_double), dimension(0:2) :: RR, E10 58 real(c_double) :: L10, D10 59 60contains !****************************************************************************************** 61 62 subroutine PairWeight1 ( W, E1_1, E1_2, E2_1, E2_2, R2_1, R2_2 ) !!!!!!!!!!!!!!!!!!!!!!!!!!! 63 real(c_double), intent(out) :: W 64 real(c_double), dimension(0:2), intent(out) :: E1_1, E1_2, E2_1, E2_2 65 real(c_double), dimension(0:2), intent(in) :: R2_1, R2_2 66 !------------------------------------------------------------------------------------------- 67 real(c_double) :: D, L20, D20, t, dWdD 68 real(c_double), dimension(0:2) :: E, E20 69 !------------------------------------------------------------------------------------------- 70 E = 0.5d+00 * ( R2_1 + R2_2 ) - RR 71 call ApplyPeriodicBC ( E ) 72 D = E(0) * E(0) + E(1) * E(1) + E(2) * E(2) 73 if ( D < D10 * D10 ) then 74 W = 1.0d+00 75 E1_1 = 0.0d+00 76 E1_2 = 0.0d+00 77 E2_1 = 0.0d+00 78 E2_2 = 0.0d+00 79 return 80 end if 81 E20 = 0.5d+00 * ( R2_2 - R2_1 ) 82 L20 = sqrt ( S_V3xx ( E20 ) + sqr ( TPMR2 ) ) 83 D20 = L10 + L20 + TPBRcutoff + RSkin 84 if ( D > D20 * D20 ) then 85 W = 0.0d+00 86 E1_1 = 0.0d+00 87 E1_2 = 0.0d+00 88 E2_1 = 0.0d+00 89 E2_2 = 0.0d+00 90 return 91 end if 92 D = sqrt ( D ) 93 E = E / D 94 E20 = E20 / L20 95 D20 = D20 - D10 96 t = ( D - D10 ) / D20 97 W = 1.0d+00 - t * t * ( 3.0d+00 - 2.0d+00 * t ) 98 dWdD = 3.0d+00 * t * ( t - 1.0d+00 ) / D20 99 E1_1 = dWdD * ( t * E10 - E ) 100 E1_2 = dWdD * ( - t * E10 - E ) 101 E2_1 = dWdD * ( E + t * E20 ) 102 E2_2 = dWdD * ( E - t * E20 ) 103 end subroutine PairWeight1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 104 105 integer(c_int) function EndWeight1 ( W, E1_1, E1_2, E2_1, E2_2, R1_1, R1_2, R2_1, R2_2 ) !!! 106 real(c_double), intent(out) :: W 107 real(c_double), dimension(0:2), intent(out) :: E1_1, E1_2, E2_1, E2_2 108 real(c_double), dimension(0:2), intent(in) :: R1_1, R1_2, R2_1, R2_2 109 !------------------------------------------------------------------------------------------- 110 real(c_double) :: D, L20 111 real(c_double) :: D1, D2, t, dWdD 112 real(c_double), dimension(0:2) :: RR, E, E20 113 !------------------------------------------------------------------------------------------- 114 E = 0.5d+00 * ( R2_1 + R2_2 - ( R1_1 + R1_2 ) ) 115 call ApplyPeriodicBC ( E ) 116 D = S_V3norm3 ( E ) 117 E20 = 0.5d+00 * ( R2_2 - R2_1 ) 118 L20 = sqrt ( S_V3xx ( E20 ) + sqr ( TPMR2 ) ) 119 D1 = L10 + L20 + TPBRcutoff + RSkin 120 if ( D < D1 ) then 121 EndWeight1 = 0 122 W = 1.0d+00 123 E1_1 = 0.0d+00 124 E1_2 = 0.0d+00 125 E2_1 = 0.0d+00 126 E2_2 = 0.0d+00 127 return 128 end if 129 D2 = D1 + TPMC3 130 if ( D > D2 ) then 131 EndWeight1 = 2 132 W = 0.0d+00 133 E1_1 = 0.0d+00 134 E1_2 = 0.0d+00 135 E2_1 = 0.0d+00 136 E2_2 = 0.0d+00 137 return 138 end if 139 EndWeight1 = 1 140 E = E / D 141 E20 = E20 / L20 142 t = ( D - D1 ) / TPMC3 143 W = 1.0d+00 - t * t * ( 3.0d+00 - 2.0d+00 * t ) 144 dWdD = 3.0d+00 * t * ( t - 1.0d+00 ) / TPMC3 145 E1_1 = dWdD * ( E10 - E ) 146 E1_2 = dWdD * ( - E10 - E ) 147 E2_1 = dWdD * ( E + E20 ) 148 E2_2 = dWdD * ( E - E20 ) 149 end function EndWeight1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 150 151 integer(c_int) function TPMInteractionFC1 ( Q, U, F1, F2, P1, P2, Pe, Pe1, R1, R2, Q1, Q2, Qe, Qe1, EType ) 152 real(c_double), intent(out) :: Q, U 153 real(c_double), dimension(0:2), intent(out) :: F1, F2, P1, P2, Pe, Pe1 154 real(c_double), dimension(0:2), intent(in) :: R1, R2, Q1, Q2, Qe, Qe1 155 integer(c_int), intent(in) :: EType 156 !------------------------------------------------------------------------------------------- 157 real(c_double), dimension(0:2) :: M, QX, Me, F1a, F2a, P1a, P2a, F1b, F2b, P1b, P2b, ER1, ER2, EQe, EQe1 158 real(c_double) :: W, W1, D, Qa, Qb, Ua, Ub, L, Pee, Peea, Peeb, DU 159 integer(c_int) :: IntSigna, IntSignb, CaseID 160 !------------------------------------------------------------------------------------------- 161 if ( EType == 0 ) then 162 TPMInteractionFC1 = TPMInteractionF ( Q, U, F1, F2, P1, P2, Pee, R1, R2, Q1, Q2, 0 ) 163 Pe = 0.0d+00 164 Pe1 = 0.0d+00 165 else if ( EType < 3 ) then 166 QX = 0.5d+00 * ( Q1 + Q2 ) 167 M = Q2 - Q1 168 L = S_V3norm3 ( M ) 169 M = M / L 170 Me = Qe - QX 171 D = S_V3norm3 ( Me ) 172 if ( EType == 1 ) then 173 TPMInteractionFC1 = TPMInteractionF ( Q, U, F1, F2, P1, P2, Pee, R1, R2, QX - D * M, QX, 1 ) 174 else 175 TPMInteractionFC1 = TPMInteractionF ( Q, U, F1, F2, P1, P2, Pee, R1, R2, QX, QX + D * M, 2 ) 176 end if 177 call TPMSegmentForces ( P1, P2, F1, F2, R1, R2, QX, M, L ) 178 Pe = ( Pee / D ) * Me 179 Pe1 = 0.0d+00 180 QX = 0.5d+00 * Pe 181 P1 = P1 + QX 182 P2 = P2 + QX 183 else 184 CaseID = EndWeight1 ( W, ER1, ER2, EQe, Eqe1, R1, R2, Qe, Qe1 ) 185 if ( CaseID < 2 ) then 186 QX = 0.5d+00 * ( Q1 + Q2 ) 187 M = Q2 - Q1 188 L = S_V3norm3 ( M ) 189 M = M / L 190 Me = Qe - QX 191 D = S_V3norm3 ( Me ) 192 if ( EType == 3 ) then 193 IntSigna = TPMInteractionF ( Qa, Ua, F1a, F2a, P1a, P2a, Peea, R1, R2, QX - D * M, QX, 1 ) 194 else 195 IntSigna = TPMInteractionF ( Qa, Ua, F1a, F2a, P1a, P2a, Peea, R1, R2, QX, QX + D * M, 2 ) 196 end if 197 call TPMSegmentForces ( P1a, P2a, F1a, F2a, R1, R2, QX, M, L ) 198 end if 199 200 if ( CaseID > 0 ) then 201 IntSignb = TPMInteractionF ( Qb, Ub, F1b, F2b, P1b, P2b, Peeb, R1, R2, Q1, Q2, 0 ) 202 end if 203 204 if ( CaseID == 0 ) then 205 TPMInteractionFC1 = IntSigna 206 Q = Qa 207 U = Ua 208 F1 = F1a 209 F2 = F2a 210 Pe = ( Peea / D ) * Me 211 Pe1 = 0.0d+00 212 QX = 0.5d+00 * Pe 213 P1 = P1a + QX 214 P2 = P2a + QX 215 else if ( CaseID == 2 ) then 216 TPMInteractionFC1 = IntSignb 217 Q = Qb 218 U = Ub 219 F1 = F1b 220 F2 = F2b 221 P1 = P1b 222 P2 = P2b 223 Pe = 0.0d+00 224 Pe1 = 0.0d+00 225 else 226 TPMInteractionFC1 = 0 227 if ( IntSigna > 0 .or. IntSignb > 0 ) TPMInteractionFC1 = 1 228 W1 = 1.0d+00 - W 229 DU = Ub - Ua 230 Q = W * Qa + W1 * Qb 231 U = W * Ua + W1 * Ub 232 Pe = ( W * Peea / D ) * Me 233 QX = 0.5d+00 * Pe 234 F1 = W * F1a + W1 * F1b + DU * ER1 235 F2 = W * F2a + W1 * F2b + DU * ER2 236 P1 = W * P1a + W1 * P1b + QX 237 P2 = W * P2a + W1 * P2b + QX 238 Pe = Pe - DU * EQe 239 Pe1 = - DU * EQe1 240 end if 241 end if 242 end function TPMInteractionFC1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 243 244 integer(c_int) function TPMInteractionFW1 ( QQ, U, U1, U2, UU, F1, F2, F, Fe, G1, G2, R1, R2, N, NMAX, R, Re, EType ) 245 real(c_double), intent(out) :: U, U1, U2 246 integer(c_int), intent(in) :: N, NMAX, EType 247 real(c_double), dimension(0:NMAX-1), intent(out) :: QQ, UU 248 real(c_double), dimension(0:2), intent(out) :: F1, F2, Fe 249 real(c_double), dimension(0:2,0:NMAX-1), intent(out) :: F, G1, G2 250 real(c_double), dimension(0:2), intent(in) :: R1, R2, Re 251 real(c_double), dimension(0:2,0:NMAX-1), intent(in) :: R 252 !------------------------------------------------------------------------------------------- 253 integer(c_int) :: i, j 254 real(c_double) :: Q, WW, DD 255 !------------------------------------------------------------------------------------------- 256 Q1 = 0.0d+00 257 Q2 = 0.0d+00 258 WW = 0.0d+00 259 Z1 = 0.0d+00 260 Z2 = 0.0d+00 261 TPMInteractionFW1 = 0 262 E10 = 0.5d+00 * ( R2 - R1 ) 263 L10 = sqrt ( S_V3xx ( E10 ) + sqr ( TPMR1 ) ) 264 D10 = TPMR1 + TPMR2 + TPMC123 * TPBRcutoff + RSkin 265 E10 = E10 / L10 266 RR = 0.5d+00 * ( R1 + R2 ) 267 do i = 0, N - 2 268 call PairWeight1 ( W(i), E1(0:2,i), E2(0:2,i), EE1(0:2,i), EE2(0:2,i), R(0:2,i), R(0:2,i+1) ) 269 Q1 = Q1 + W(i) * R(0:2,i) 270 Q2 = Q2 + W(i) * R(0:2,i+1) 271 WW = WW + W(i) 272 Z1 = Z1 + E1(0:2,i) 273 Z2 = Z2 + E2(0:2,i) 274 end do 275 if ( WW .le. TPGeomPrec ) return 276 Q1 = Q1 / WW 277 Q2 = Q2 / WW 278 Z1 = Z1 / WW 279 Z2 = Z2 / WW 280 if ( EType == 1 ) then 281 Qe = R(0:2,0) 282 Qe1 = R(0:2,1) 283 else if ( EType == 2 ) then 284 Qe = R(0:2,N-1) 285 Qe1 = R(0:2,N-2) 286 else if ( EType == 3 ) then 287 Qe = Re 288 Qe1 = R(0:2,0) 289 else if ( EType == 4 ) then 290 Qe = Re 291 Qe1 = R(0:2,N-1) 292 else 293 Qe = 0.0d+00 294 Qe1 = 0.0d+00 295 end if 296 297 TPMInteractionFW1 = TPMInteractionFC1 ( Q, U, F1, F2, S1, S2, Pe, Pe1, R1, R2, Q1, Q2, Qe, Qe1, EType ) 298 if ( TPMInteractionFW1 == 0 ) return 299 300 W(0:N-2) = W(0:N-2) / WW 301 E1(0:2,0:N-2) = E1(0:2,0:N-2) / WW 302 E2(0:2,0:N-2) = E2(0:2,0:N-2) / WW 303 EE1(0:2,0:N-2) = EE1(0:2,0:N-2) / WW 304 EE2(0:2,0:N-2) = EE2(0:2,0:N-2) / WW 305 G1(0:2,0:N-1) = 0.0d+00 306 G2(0:2,0:N-1) = 0.0d+00 307 U1 = 0.25d+00 * U 308 U2 = U1 309 UU = 0.0d+00 310 do i = 0, N - 2 311 QQ(i) = W(i) * Q 312 DD = W(i) * U1 313 UU(i) = UU(i) + DD 314 UU(i+1) = UU(i+1) + DD 315 end do 316 do i = 0, N - 2 317 C(i) = S_V3xV3 ( S1, R(0:2,i) ) + S_V3xV3 ( S2, R(0:2,i+1) ) 318 F1 = F1 + C(i) * ( E1(0:2,i) - W(i) * Z1 ) 319 F2 = F2 + C(i) * ( E2(0:2,i) - W(i) * Z2 ) 320 end do 321 F(0:2,0) = W(0) * S1 322 do j = 0, N - 2 323 if ( j == 0 ) then 324 DR = EE1(0:2,0) * ( 1.0d+00 - W(0) ) 325 else 326 DR = - W(j) * EE1(0:2,0) 327 end if 328 F(0:2,0) = F(0:2,0) + C(j) * DR 329 end do 330 do i = 1, N - 2 331 G1(0:2,i) = W(i-1) * S2 332 G2(0:2,i) = W(i) * S1 333 do j = 0, N - 2 334 if ( j == i ) then 335 G1(0:2,i) = G1(0:2,i) - C(j) * W(j) * EE2(0:2,i-1) 336 G2(0:2,i) = G2(0:2,i) + C(j) * ( EE1(0:2,j) - W(j) * EE1(0:2,i) ) 337 else if ( j == i - 1 ) then 338 G1(0:2,i) = G1(0:2,i) + C(j) * ( EE2(0:2,j) - W(j) * EE2(0:2,i-1) ) 339 G2(0:2,i) = G2(0:2,i) - C(j) * W(j) * EE1(0:2,i) 340 else 341 G1(0:2,i) = G1(0:2,i) - C(j) * W(j) * EE2(0:2,i-1) 342 G2(0:2,i) = G2(0:2,i) - C(j) * W(j) * EE1(0:2,i) 343 end if 344 end do 345 F(0:2,i) = G1(0:2,i) + G2(0:2,i) 346 end do 347 F(0:2,N-1) = W(N-2) * S2 348 do j = 0, N - 2 349 if ( j == N - 2 ) then 350 DR = EE2(0:2,N-2) * ( 1.0d+00 - W(N-2) ) 351 else 352 DR = - W(j) * EE2(0:2,N-2) 353 end if 354 F(0:2,N-1) = F(0:2,N-1) + C(j) * DR 355 end do 356 Fe = 0.0d+00 357 if ( EType == 1 ) then 358 F(0:2,0) = F(0:2,0) - Pe 359 else if ( EType == 2 ) then 360 F(0:2,N-1) = F(0:2,N-1) - Pe 361 else if ( EType == 3 ) then 362 F(0:2,0) = F(0:2,0) - Pe1 363 Fe = - Pe 364 else if ( EType == 4 ) then 365 F(0:2,N-1) = F(0:2,N-1) - Pe1 366 Fe = - Pe 367 end if 368 G1(0:2,N-1) = F(0:2,N-1) 369 G2(0:2,0) = F(0:2,0) 370 end function TPMInteractionFW1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 371 372end module TPMM1 !********************************************************************************** 373