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