ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/gb.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/gb.F90 (file contents):
Revision 2204 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 2367 by kdaily, Fri Oct 14 15:44:37 2005 UTC

# Line 44 | Line 44 | module gb_pair
44    use force_globals
45    use definitions
46    use simulation
47 +  use atype_module
48 +  use vector_class
49 +  use status
50 +  use lj
51   #ifdef IS_MPI
52    use mpiSimulation
53   #endif
54 <
54 >  
55    implicit none
56  
57 <  PRIVATE
57 >  private
58  
59 +  logical, save :: haveGayBerneLJMap = .false.
60    logical, save :: gb_pair_initialized = .false.
61 +  logical, save :: gb_lj_pair_initialized = .false.
62    real(kind=dp), save :: gb_sigma
63    real(kind=dp), save :: gb_l2b_ratio
64    real(kind=dp), save :: gb_eps
65    real(kind=dp), save :: gb_eps_ratio
66    real(kind=dp), save :: gb_mu
67    real(kind=dp), save :: gb_nu
68 +  real(kind=dp), save :: lj_sigma
69 +  real(kind=dp), save :: lj_eps
70 +  real(kind=dp), save :: gb_sigma_l
71 +  real(kind=dp), save :: gb_eps_l
72  
73    public :: check_gb_pair_FF
74    public :: set_gb_pair_params
75    public :: do_gb_pair
76 +  public :: getGayBerneCut
77 + !!$  public :: check_gb_lj_pair_FF
78 + !!$  public :: set_gb_lj_pair_params
79 +  public :: do_gb_lj_pair
80 +  public :: createGayBerneLJMap
81 +  public :: destroyGayBerneTypes
82 +  public :: complete_GayBerne_FF
83 + !!may not need
84 +  type, private :: LJtype
85 +     integer  :: atid
86 +     real(kind=dp) :: sigma
87 +     real(kind=dp) :: epsilon
88 +  end type LJtype
89 + !!may not need
90 +  type, private :: LJList
91 +     integer       :: Nljtypes =0
92 +     integer       :: currentLJtype= 0
93 +     type(LJtype), pointer :: LJtype(:) =>  null()
94 +     integer, pointer      :: atidToLJtype(:) =>null()
95 +  end type LJList
96  
97 +  type(LJList), save :: LJMap
98 +  
99 +  type :: GayBerneLJ
100 + !!$     integer :: atid
101 + !!$     real(kind = dp ),pointer, dimension(:) ::  epsil_GB      =>null()
102 + !!$     real(kind = dp ),pointer, dimension(:) ::  sigma_GB        =>null()
103 + !!$     real(kind = dp ),pointer, dimension(:) ::  epsilon_ratio   =>null()
104 + !!$     real(kind = dp ),pointer, dimension(:) ::  sigma_ratio     =>null()
105 + !!$     real(kind = dp ),pointer, dimension(:) ::  mu              =>null()
106 +    
107 +     real(kind = dp ) :: sigma_l
108 +     real(kind = dp ) :: sigma_s
109 +     real(kind = dp ) :: sigma_ratio
110 +     real(kind = dp ) :: eps_s
111 +     real(kind = dp ) :: eps_l
112 +     real(kind = dp ) :: eps_ratio
113 +     integer          :: c_ident
114 +     integer          :: status
115 +  end type GayBerneLJ
116 +
117 + !!$  type, private :: gayberneLJlist
118 + !!$     integer:: n_gaybernelj = 0
119 + !!$     integer:: currentgayberneLJ = 0
120 + !!$     type(GayBerneLJ),pointer :: GayBerneLJ(:) => null()
121 + !!$     integer, pointer         :: atidToGayBerneLJ(:) => null()
122 + !!$  end type gayberneLJlist
123 +
124 +  type(gayberneLJ), dimension(:), allocatable :: gayBerneLJMap
125 +
126   contains
127  
128    subroutine check_gb_pair_FF(status)
# Line 73 | Line 132 | contains
132      return
133    end subroutine check_gb_pair_FF
134  
135 + !!$  subroutine check_gb_lj_pair_FF(status)
136 + !!$    integer :: status
137 + !!$    status = -1
138 + !!$    if (gb_lj_pair_initialized) status = 0
139 + !!$    return
140 + !!$  end subroutine check_gb_lj_pair_FF
141 +
142    subroutine set_gb_pair_params(sigma, l2b_ratio, eps, eps_ratio, mu, nu)
143      real( kind = dp ), intent(in) :: sigma, l2b_ratio, eps, eps_ratio
144      real( kind = dp ), intent(in) :: mu, nu
145 <
145 >    integer :: ntypes, nljtypes
146 > !!    integer, intent(in) :: c_ident
147 >    integer, pointer :: MatchList(:) => null ()
148 >    integer :: status
149      gb_sigma = sigma
150      gb_l2b_ratio = l2b_ratio
151      gb_eps = eps
152      gb_eps_ratio = eps_ratio
153      gb_mu = mu
154      gb_nu = nu
155 +    gb_sigma_l = gb_sigma*l2b_ratio
156 +    gb_eps_l = gb_eps*gb_eps_ratio
157 +    status = 0
158  
159 <    gb_pair_initialized = .true.
159 >
160 >      
161 >
162      return
163    end subroutine set_gb_pair_params
164 +  
165 + !!$  subroutine set_gb_lj_pair_params(sigma_gb, l2b_ratio, eps_gb, eps_ratio, mu, nu, sigma_lj, eps_lj, c_ident, status)
166 + !!$    real( kind = dp ), intent(in) :: sigma_gb, l2b_ratio, eps_gb, eps_ratio, mu, nu
167 + !!$    real( kind = dp ), intent(in) :: sigma_lj, eps_lj
168 + !!$    integer, intent(in) :: c_ident
169 + !!$    integer :: nLJTYPES, nGayBerneTypes, ntypes, current, status
170 + !!$
171 + !!$    integer :: myATID
172 + !!$    logical :: thisProperty
173 + !!$    real(kind=dp):: fake
174 + !!$    
175 + !!$    status = 0
176 + !!$
177 + !!$    if(.not.associated(LJMap%Ljtype)) then
178 +
179 + !!$       call getMatchingElementList(atypes, "is_GayBerne", .true., &
180 + !!$            nGayBerneTypes, MatchList)
181 + !!$      
182 + !!$       call getMatchingElementList(atypes, "is_LennardJones", .true., &
183 + !!$            nLJTypes, MatchList)
184 + !!$
185 + !!$       LJMap%nLJtypes = nLJTypes
186 + !!$
187 + !!$       allocate(LJMap% LJtype(nLJTypes))
188 + !!$
189 + !!$       ntypes = getSize(atypes)
190 + !!$
191 + !!$       allocate(LJMap%atidToLJtype(ntypes))
192 + !!$
193 + !!$   endif
194 + !!$
195 + !!$   LJmap%currentLJtype =  LJMap%currentLJtype + 1
196 + !!$
197 + !!$   current = LJMap%currentLJtype
198 + !!$   LJMap%atidToLJtype(myATID)    = current
199 + !!$   myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
200 + !!$   call getElementProperty(atypes, myATID, "is_LennardJones",thisProperty)
201 + !!$   if(thisProperty) then
202 + !!$
203 + !!$      LJMap%LJtype(current)%atid    = myatid
204 + !!$!!for testing
205 + !!$      fake = getSigma(myATID)
206 + !!$      LJMap%LJtype(current)%sigma   = getSigma(myATID)
207 + !!$      LJMap%LJtype(current)%epsilon = getEpsilon(myATID)
208 + !!$   end if
209 +
210 + !!$    gb_sigma = sigma_gb
211 + !!$    gb_l2b_ratio = l2b_ratio
212 + !!$    gb_eps = eps_gb
213 + !!$    gb_eps_ratio = eps_ratio
214 + !!    gb_mu = mu
215 + !!    gb_nu = nu
216 + !!$    
217 + !!$    lj_sigma = sigma_lj
218 + !!$    lj_eps = eps_lj
219 +
220 + !!    gb_lj_pair_initialized = .true.
221 + !!$  endsubroutine set_gb_lj_pair_params
222 +
223 +  subroutine createGayBerneLJMap
224 +    integer :: ntypes, i, j
225 +    real(kind=dp) :: s1, s2, e1, e2
226 +    real(kind=dp) :: sigma_s,sigma_l,eps_s, eps_l
227 +    
228 +    if(LJMap%currentLJtype == 0)then
229 +       call handleError("gayberneLJ", "no members in gayberneLJMap")
230 +       return
231 +    end if
232 +
233 +    ntypes = getSize(atypes)
234 +    
235 +    if(.not.allocated(gayBerneLJMap))then
236 +       allocate(gayBerneLJMap(ntypes))
237 +    endif
238 +    
239 +    do i = 1, ntypes
240 +       s1 = LJMap%LJtype(i)%sigma
241 +       e1 = LJMap%LJtype(i)%epsilon
242 +    
243 + !!$       sigma_s = 0.5d0 *(s1+gb_sigma)
244 + !!$       sigma_l = 0.5d0 *(s1+gb_sigma*gb_l2b_ratio)
245 +       sigma_s = gb_sigma
246 +       sigma_l = gb_sigma*gb_l2b_ratio
247 +       gayBerneLJMap(i)%sigma_s = sigma_s
248 +       gayBerneLJMap(i)%sigma_l = sigma_l
249 +       gayBerneLJMap(i)%sigma_ratio = sigma_l/sigma_s
250 +       eps_s = dsqrt(e1*gb_eps)
251 +       eps_l = dsqrt(e1*gb_eps_l)
252 +       gayBerneLJMap(i)%eps_s = eps_s
253 +       gayBerneLJMap(i)%eps_l = eps_l
254 +       gayBerneLJMap(i)%eps_ratio = eps_l/eps_s
255 +    enddo
256 +    haveGayBerneLJMap = .true.
257 +    gb_lj_pair_initialized = .true.
258 +  endsubroutine createGayBerneLJMap
259 +  !! gay berne cutoff should be a parameter in globals, this is a temporary
260 +  !! work around - this should be fixed when gay berne is up and running
261 +  function getGayBerneCut(atomID) result(cutValue)
262 +    integer, intent(in) :: atomID !! nah... we don't need to use this...
263 +    real(kind=dp) :: cutValue
264  
265 +    cutValue = gb_l2b_ratio*gb_sigma*2.5_dp
266 +  end function getGayBerneCut
267  
268    subroutine do_gb_pair(atom1, atom2, d, r, r2, sw, vpair, fpair, &
269         pot, A, f, t, do_pot)
270 <
270 >    
271      integer, intent(in) :: atom1, atom2
272      integer :: id1, id2
273      real (kind=dp), intent(inout) :: r, r2
# Line 132 | Line 308 | contains
308      real(kind=dp) :: term2a, term2b, term2u1x, term2u1y, term2u1z
309      real(kind=dp) :: term2u2x, term2u2y, term2u2z
310      real(kind=dp) :: yick1, yick2, mess1, mess2
311 <
311 >    
312      s2 = (gb_l2b_ratio)**2
313      emu = (gb_eps_ratio)**(1.0d0/gb_mu)
314  
# Line 157 | Line 333 | contains
333      ul2(1) = A(3,atom2)
334      ul2(2) = A(6,atom2)
335      ul2(3) = A(9,atom2)
160 – #endif
336  
337 + #endif
338 +    
339      dru1dx = ul1(1)
340      dru2dx = ul2(1)
341      dru1dy = ul1(2)
342      dru2dy = ul2(2)
343      dru1dz = ul1(3)
344      dru2dz = ul2(3)
345 +        
346  
347 +
348      drdx = d(1) / r
349      drdy = d(2) / r
350      drdz = d(3) / r
351 <
351 >    
352      ! do some dot products:
353      ! NB the r in these dot products is the actual intermolecular vector,
354      ! and is not the unit vector in that direction.
355 <
355 >    
356      rdotu1 = d(1)*ul1(1) + d(2)*ul1(2) + d(3)*ul1(3)
357      rdotu2 = d(1)*ul2(1) + d(2)*ul2(2) + d(3)*ul2(3)
358      u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) +  ul1(3)*ul2(3)
# Line 184 | Line 363 | contains
363      ! We note however, that there are some major typos in that Appendix
364      ! of the Luckhurst paper, particularly in equations A23, A29 and A31
365      ! We have attempted to correct them below.
366 <
366 >    
367      dotsum = rdotu1+rdotu2
368      dotdiff = rdotu1-rdotu2
369      ds2 = dotsum*dotsum
370      dd2 = dotdiff*dotdiff
371 <
371 >  
372      opXdot = 1.0d0 + Chi*u1dotu2
373      omXdot = 1.0d0 - Chi*u1dotu2
374      opXpdot = 1.0d0 + ChiPrime*u1dotu2
375      omXpdot = 1.0d0 - ChiPrime*u1dotu2
376 <
376 >  
377      line1a = dotsum/opXdot
378      line1bx = dru1dx + dru2dx
379      line1by = dru1dy + dru2dy
380      line1bz = dru1dz + dru2dz
381 <
381 >    
382      line2a = dotdiff/omXdot
383      line2bx = dru1dx - dru2dx
384      line2by = dru1dy - dru2dy
385      line2bz = dru1dz - dru2dz
386 <
386 >    
387      term1x = -Chi*(line1a*line1bx + line2a*line2bx)/r2
388      term1y = -Chi*(line1a*line1by + line2a*line2by)/r2
389      term1z = -Chi*(line1a*line1bz + line2a*line2bz)/r2
390 <
390 >    
391      line3a = ds2/opXdot
392      line3b = dd2/omXdot
393      line3 = Chi*(line3a + line3b)/r4
394      line3x = d(1)*line3
395      line3y = d(2)*line3
396      line3z = d(3)*line3
397 <
397 >    
398      dgdx = term1x + line3x
399      dgdy = term1y + line3y
400      dgdz = term1z + line3z
# Line 226 | Line 405 | contains
405      term1u2x = 2.0d0*(line1a-line2a)*d(1)
406      term1u2y = 2.0d0*(line1a-line2a)*d(2)
407      term1u2z = 2.0d0*(line1a-line2a)*d(3)
408 <
408 >    
409      term2a = -line3a/opXdot
410      term2b =  line3b/omXdot
411 <
411 >    
412      term2u1x = Chi*ul2(1)*(term2a + term2b)
413      term2u1y = Chi*ul2(2)*(term2a + term2b)
414      term2u1z = Chi*ul2(3)*(term2a + term2b)
415      term2u2x = Chi*ul1(1)*(term2a + term2b)
416      term2u2y = Chi*ul1(2)*(term2a + term2b)
417      term2u2z = Chi*ul1(3)*(term2a + term2b)
418 <
418 >    
419      pref = -Chi*0.5d0/r2
420  
421      dgdu1x = pref*(term1u1x+term2u1x)
# Line 247 | Line 426 | contains
426      dgdu2z = pref*(term1u2z+term2u2z)
427  
428      g = 1.0d0 - Chi*(line3a + line3b)/(2.0d0*r2)
429 <
429 >  
430      BigR = (r - gb_sigma*(g**(-0.5d0)) + gb_sigma)/gb_sigma
431      Ri = 1.0d0/BigR
432      Ri2 = Ri*Ri
# Line 261 | Line 440 | contains
440      dBigRdx = drdx/gb_sigma + dgdx*gfact
441      dBigRdy = drdy/gb_sigma + dgdy*gfact
442      dBigRdz = drdz/gb_sigma + dgdz*gfact
443 +
444      dBigRdu1x = dgdu1x*gfact
445      dBigRdu1y = dgdu1y*gfact
446      dBigRdu1z = dgdu1z*gfact
# Line 281 | Line 461 | contains
461      line3x = d(1)*line3
462      line3y = d(2)*line3
463      line3z = d(3)*line3
464 <
464 >    
465      dgpdx = term1x + line3x
466      dgpdy = term1y + line3y
467      dgpdz = term1z + line3z
468 <
469 <    term1u1x = 2.0d0*(line1a+line2a)*d(1)
470 <    term1u1y = 2.0d0*(line1a+line2a)*d(2)
471 <    term1u1z = 2.0d0*(line1a+line2a)*d(3)
468 >    
469 >    term1u1x = 2.00d0*(line1a+line2a)*d(1)
470 >    term1u1y = 2.00d0*(line1a+line2a)*d(2)
471 >    term1u1z = 2.00d0*(line1a+line2a)*d(3)
472      term1u2x = 2.0d0*(line1a-line2a)*d(1)
473      term1u2y = 2.0d0*(line1a-line2a)*d(2)
474      term1u2z = 2.0d0*(line1a-line2a)*d(3)
475  
476      term2a = -line3a/opXpdot
477      term2b =  line3b/omXpdot
478 <
478 >    
479      term2u1x = ChiPrime*ul2(1)*(term2a + term2b)
480      term2u1y = ChiPrime*ul2(2)*(term2a + term2b)
481      term2u1z = ChiPrime*ul2(3)*(term2a + term2b)
482      term2u2x = ChiPrime*ul1(1)*(term2a + term2b)
483      term2u2y = ChiPrime*ul1(2)*(term2a + term2b)
484      term2u2z = ChiPrime*ul1(3)*(term2a + term2b)
485 <
485 >  
486      pref = -ChiPrime*0.5d0/r2
487 <
487 >    
488      dgpdu1x = pref*(term1u1x+term2u1x)
489      dgpdu1y = pref*(term1u1y+term2u1y)
490      dgpdu1z = pref*(term1u1z+term2u1z)
491      dgpdu2x = pref*(term1u2x+term2u2x)
492      dgpdu2y = pref*(term1u2y+term2u2y)
493      dgpdu2z = pref*(term1u2z+term2u2z)
494 <
494 >    
495      gp = 1.0d0 - ChiPrime*(line3a + line3b)/(2.0d0*r2)
496      gmu = gp**gb_mu
497      gpi = 1.0d0 / gp
498      gmum = gmu*gpi
499  
320 –    ! write(*,*) atom1, atom2, Chi, u1dotu2
500      curlyE = 1.0d0/dsqrt(1.0d0 - Chi*Chi*u1dotu2*u1dotu2)
501 <
501 > !!$
502 > !!$    dcE = -(curlyE**3)*Chi*Chi*u1dotu2
503      dcE = (curlyE**3)*Chi*Chi*u1dotu2
504  
505      dcEdu1x = dcE*ul2(1)
# Line 328 | Line 508 | contains
508      dcEdu2x = dcE*ul1(1)
509      dcEdu2y = dcE*ul1(2)
510      dcEdu2z = dcE*ul1(3)
511 <
511 >    
512      enu = curlyE**gb_nu
513      enum = enu/curlyE
514 <
514 >  
515      eps = gb_eps*enu*gmu
516  
517      yick1 = gb_eps*enu*gb_mu*gmum
# Line 343 | Line 523 | contains
523      depsdu2x = yick1*dgpdu2x + yick2*dcEdu2x
524      depsdu2y = yick1*dgpdu2y + yick2*dcEdu2y
525      depsdu2z = yick1*dgpdu2z + yick2*dcEdu2z
526 <
526 >    
527      R126 = Ri12 - Ri6
528      R137 = 6.0d0*Ri7 - 12.0d0*Ri13
529 <
529 >    
530      mess1 = gmu*R137
531      mess2 = R126*gb_mu*gmum
532 <
532 >    
533      dUdx = 4.0d0*gb_eps*enu*(mess1*dBigRdx + mess2*dgpdx)*sw
534      dUdy = 4.0d0*gb_eps*enu*(mess1*dBigRdy + mess2*dgpdy)*sw
535      dUdz = 4.0d0*gb_eps*enu*(mess1*dBigRdz + mess2*dgpdz)*sw
536 <
536 >    
537      dUdu1x = 4.0d0*(R126*depsdu1x + eps*R137*dBigRdu1x)*sw
538      dUdu1y = 4.0d0*(R126*depsdu1y + eps*R137*dBigRdu1y)*sw
539      dUdu1z = 4.0d0*(R126*depsdu1z + eps*R137*dBigRdu1z)*sw
540      dUdu2x = 4.0d0*(R126*depsdu2x + eps*R137*dBigRdu2x)*sw
541      dUdu2y = 4.0d0*(R126*depsdu2y + eps*R137*dBigRdu2y)*sw
542      dUdu2z = 4.0d0*(R126*depsdu2z + eps*R137*dBigRdu2z)*sw
543 <
543 >      
544   #ifdef IS_MPI
545      f_Row(1,atom1) = f_Row(1,atom1) + dUdx
546      f_Row(2,atom1) = f_Row(2,atom1) + dUdy
547      f_Row(3,atom1) = f_Row(3,atom1) + dUdz
548 <
548 >    
549      f_Col(1,atom2) = f_Col(1,atom2) - dUdx
550      f_Col(2,atom2) = f_Col(2,atom2) - dUdy
551      f_Col(3,atom2) = f_Col(3,atom2) - dUdz
552 <
553 <    t_Row(1,atom1) = t_Row(1,atom1) - ul1(2)*dUdu1z + ul1(3)*dUdu1y
554 <    t_Row(2,atom1) = t_Row(2,atom1) - ul1(3)*dUdu1x + ul1(1)*dUdu1z
555 <    t_Row(3,atom1) = t_Row(3,atom1) - ul1(1)*dUdu1y + ul1(2)*dUdu1x
556 <
557 <    t_Col(1,atom2) = t_Col(1,atom2) - ul2(2)*dUdu2z + ul2(3)*dUdu2y
558 <    t_Col(2,atom2) = t_Col(2,atom2) - ul2(3)*dUdu2x + ul2(1)*dUdu2z
559 <    t_Col(3,atom2) = t_Col(3,atom2) - ul2(1)*dUdu2y + ul2(2)*dUdu2x
552 >    
553 >    t_Row(1,atom1) = t_Row(1,atom1)- ul1(3)*dUdu1y + ul1(2)*dUdu1z
554 >    t_Row(2,atom1) = t_Row(2,atom1)- ul1(1)*dUdu1z + ul1(3)*dUdu1x
555 >    t_Row(3,atom1) = t_Row(3,atom1)- ul1(2)*dUdu1x + ul1(1)*dUdu1y
556 >    
557 >    t_Col(1,atom2) = t_Col(1,atom2) - ul2(3)*dUdu2y + ul2(2)*dUdu2z
558 >    t_Col(2,atom2) = t_Col(2,atom2) - ul2(1)*dUdu2z + ul2(3)*dUdu2x
559 >    t_Col(3,atom2) = t_Col(3,atom2) - ul2(2)*dUdu2x + ul2(1)*dUdu2y
560   #else
561      f(1,atom1) = f(1,atom1) + dUdx
562      f(2,atom1) = f(2,atom1) + dUdy
563      f(3,atom1) = f(3,atom1) + dUdz
564 <
564 >    
565      f(1,atom2) = f(1,atom2) - dUdx
566      f(2,atom2) = f(2,atom2) - dUdy
567      f(3,atom2) = f(3,atom2) - dUdz
568 <
569 <    t(1,atom1) = t(1,atom1) - ul1(2)*dUdu1z + ul1(3)*dUdu1y
570 <    t(2,atom1) = t(2,atom1) - ul1(3)*dUdu1x + ul1(1)*dUdu1z
571 <    t(3,atom1) = t(3,atom1) - ul1(1)*dUdu1y + ul1(2)*dUdu1x
572 <
573 <    t(1,atom2) = t(1,atom2) - ul2(2)*dUdu2z + ul2(3)*dUdu2y
574 <    t(2,atom2) = t(2,atom2) - ul2(3)*dUdu2x + ul2(1)*dUdu2z
575 <    t(3,atom2) = t(3,atom2) - ul2(1)*dUdu2y + ul2(2)*dUdu2x
568 >    
569 >    t(1,atom1) = t(1,atom1)- ul1(3)*dUdu1y + ul1(2)*dUdu1z
570 >    t(2,atom1) = t(2,atom1)- ul1(1)*dUdu1z + ul1(3)*dUdu1x
571 >    t(3,atom1) = t(3,atom1)- ul1(2)*dUdu1x + ul1(1)*dUdu1y
572 >    
573 >    t(1,atom2) = t(1,atom2)- ul2(3)*dUdu2y + ul2(2)*dUdu2z
574 >    t(2,atom2) = t(2,atom2)- ul2(1)*dUdu2z + ul2(3)*dUdu2x
575 >    t(3,atom2) = t(3,atom2)- ul2(2)*dUdu2x + ul2(1)*dUdu2y
576   #endif
577 <
577 >  
578      if (do_pot) then
579   #ifdef IS_MPI
580         pot_row(atom1) = pot_row(atom1) + 2.0d0*eps*R126*sw
# Line 412 | Line 592 | contains
592      id1 = atom1
593      id2 = atom2
594   #endif
595 <
595 >    
596      if (molMembershipList(id1) .ne. molMembershipList(id2)) then
597 <
597 >      
598         fpair(1) = fpair(1) + dUdx
599         fpair(2) = fpair(2) + dUdy
600         fpair(3) = fpair(3) + dUdz
601 +      
602 +    endif
603 +    
604 +    return
605 +  end subroutine do_gb_pair
606  
607 +  subroutine do_gb_lj_pair(atom1, atom2, d, r, r2, sw, vpair, fpair, &
608 +       pot, A, f, t, do_pot)
609 +    
610 +    integer, intent(in) :: atom1, atom2
611 +    integer :: id1, id2
612 +    real (kind=dp), intent(inout) :: r, r2
613 +    real (kind=dp), dimension(3), intent(in) :: d
614 +    real (kind=dp), dimension(3), intent(inout) :: fpair
615 +    real (kind=dp) :: pot, sw, vpair
616 +    real (kind=dp), dimension(9,nLocal) :: A
617 +    real (kind=dp), dimension(3,nLocal) :: f
618 +    real (kind=dp), dimension(3,nLocal) :: t
619 +    logical, intent(in) :: do_pot
620 +    real (kind = dp), dimension(3) :: ul
621 +
622 + !!  real(kind=dp) :: lj2, s_lj2pperp2,s_perpppar2,eabnu, epspar
623 +    real(kind=dp) :: spar, sperp, slj, par2, perp2, sc, slj2
624 +    real(kind=dp) :: s_par2mperp2, s_lj2ppar2
625 +    real(kind=dp) :: enot, eperp, epar, eab, eabf,moom, mum1
626 + !!    real(kind=dp) :: e_ljpperp, e_perpmpar, e_ljppar
627 +    real(kind=dp) :: dx, dy, dz, drdx,drdy,drdz, rdotu
628 + !!    real(kind=dp) :: ct, dctdx, dctdy, dctdz, dctdux, dctduy, dctduz
629 +    real(kind=dp) :: mess, sab, dsabdct, eprime, deprimedct, depmudct
630 +    real(kind=dp) :: epmu, depmudx, depmudy, depmudz
631 +    real(kind=dp) :: depmudux, depmuduy, depmuduz
632 +    real(kind=dp) :: BigR, dBigRdx, dBigRdy, dBigRdz
633 +    real(kind=dp) :: dBigRdux, dBigRduy, dBigRduz
634 +    real(kind=dp) :: dUdx, dUdy, dUdz, dUdux, dUduy, dUduz, e0
635 +    real(kind=dp) :: Ri, Ri3, Ri6, Ri7, Ril2, Ri13, Rl26, R137, prefactor
636 +    real(kind=dp) :: chipoalphap2, chioalpha2, ec, epsnot
637 +    real(kind=dp) :: drdotudx, drdotudy, drdotudz    
638 +    real(kind=dp) :: ljeps, ljsigma, sigmaratio, sigmaratioi
639 +    integer :: ljt1, ljt2, atid1, atid2
640 +    logical :: thisProperty
641 + #ifdef IS_MPI
642 +    atid1 = atid_Row(atom1)
643 +    atid2 = atid_Col(atom2)
644 + #else
645 +    atid1 = atid(atom1)
646 +    atid2 = atid(atom2)
647 + #endif
648 +    ri =1/r
649 +    
650 +    dx = d(1)
651 +    dy = d(2)
652 +    dz = d(3)
653 +    
654 +    drdx = dx *ri
655 +    drdy = dy *ri
656 +    drdz = dz *ri
657 +  
658 +  
659 +
660 +    if(.not.haveGayBerneLJMap) then
661 +       call createGayBerneLJMap()
662      endif
663 + !!$   write(*,*) "in gbljpair"
664 +    call getElementProperty(atypes, atid1, "is_LennardJones",thisProperty)
665 + !!$    write(*,*) thisProperty
666 +    if(thisProperty.eqv..true.)then
667 + #ifdef IS_MPI
668 +       ul(1) = A_Row(3,atom2)
669 +       ul(2) = A_Row(6,atom2)
670 +       ul(3) = A_Row(9,atom2)
671  
672 + #else
673 +       ul(1) = A(3,atom2)
674 +       ul(2) = A(6,atom2)
675 +       ul(3) = A(9,atom2)
676 + #endif
677 +
678 +       rdotu = (dx*ul(1)+dy*ul(2)+dz*ul(3))*ri
679 +
680 +       ljt1 = LJMap%atidtoLJtype(atid1)
681 +       ljt2 = LJMap%atidtoLJtype(atid2)
682 +      
683 +       ljeps =  LJMap%LJtype(ljt1)%epsilon
684 + !!$       write(*,*) "ljeps"
685 + !!$       write(*,*) ljeps
686 +       drdotudx = ul(1)*ri-rdotu*dx*ri*ri
687 +       drdotudy = ul(2)*ri-rdotu*dy*ri*ri
688 +       drdotudz = ul(3)*ri-rdotu*dz*ri*ri
689 +
690 +    
691 +       moom =  1.0d0 / gb_mu
692 +       mum1 = gb_mu-1
693 +
694 +       sperp = GayBerneLJMap(ljt1)%sigma_s
695 +       spar =  GayBerneLJMap(ljt1)%sigma_l
696 +       slj = LJMap%LJtype(ljt1)%sigma
697 +       slj2 = slj*slj
698 + !!$       write(*,*) "spar"
699 + !!$       write(*,*) sperp
700 + !!$       write(*,*) spar
701 + !!    chioalpha2 = s_par2mperp2/s_lj2ppar2
702 +       chioalpha2 =1-((sperp+slj)*(sperp+slj))/((spar+slj)*(spar+slj))
703 +       sc = (sperp+slj)/2.0d0
704 +  
705 +       par2 = spar*spar
706 +       perp2 = sperp*sperp
707 +       s_par2mperp2 = par2 - perp2
708 +       s_lj2ppar2 = slj2 + par2
709 +
710 +
711 + !! check these ! left from old code
712 + !! kdaily e0 may need to be (gb_eps + lj_eps)/2
713 +    
714 +       eperp = dsqrt(gb_eps*ljeps)
715 +       epar = eperp*gb_eps_ratio
716 +       enot = dsqrt(ljeps*eperp)
717 +       chipoalphap2 = 1+(dsqrt(epar*ljeps)/enot)**moom
718 + !! to so mess matchs cleaver (eq 20)
719 +    
720 +       mess = 1-rdotu*rdotu*chioalpha2
721 +       sab = 1.0d0/dsqrt(mess)
722 +       dsabdct = sc*sab*sab*sab*rdotu*chioalpha2
723 +      
724 +       eab = 1-chipoalphap2*rdotu*rdotu
725 +       eabf = enot*eab**gb_mu
726 +       depmudct = -2*enot*chipoalphap2*gb_mu*rdotu*eab**mum1
727 +      
728 +      
729 +       BigR = (r - sab*sc + sc)/sc
730 +       dBigRdx = (drdx -dsabdct*drdotudx)/sc
731 +       dBigRdy = (drdy -dsabdct*drdotudy)/sc
732 +       dBigRdz = (drdz -dsabdct*drdotudz)/sc
733 +       dBigRdux = (-dsabdct*drdx)/sc
734 +       dBigRduy = (-dsabdct*drdy)/sc
735 +       dBigRduz = (-dsabdct*drdz)/sc
736 +      
737 +       depmudx = depmudct*drdotudx
738 +       depmudy = depmudct*drdotudy
739 +       depmudz = depmudct*drdotudz
740 +       depmudux = depmudct*drdx
741 +       depmuduy = depmudct*drdy
742 +       depmuduz = depmudct*drdz
743 +      
744 +       Ri = 1.0d0/BigR
745 +       Ri3 = Ri*Ri*Ri
746 +       Ri6 = Ri3*Ri3
747 +       Ri7 = Ri6*Ri
748 +       Ril2 = Ri6*Ri6
749 +       Ri13 = Ri6*Ri7
750 +       Rl26 = Ril2 - Ri6
751 +       R137 = 6.0d0*Ri7 - 12.0d0*Ri13
752 +      
753 +       prefactor = 4.0d0
754 +      
755 +       dUdx = prefactor*(eabf*R137*dBigRdx + Rl26*depmudx)
756 +       dUdy = prefactor*(eabf*R137*dBigRdy + Rl26*depmudy)
757 +       dUdz = prefactor*(eabf*R137*dBigRdz + Rl26*depmudz)
758 +       dUdux = prefactor*(eabf*R137*dBigRdux + Rl26*depmudux)
759 +       dUduy = prefactor*(eabf*R137*dBigRduy + Rl26*depmuduy)
760 +       dUduz = prefactor*(eabf*R137*dBigRduz + Rl26*depmuduz)
761 +      
762 + #ifdef IS_MPI
763 +       f_Row(1,atom1) = f_Row(1,atom1) - dUdx
764 +       f_Row(2,atom1) = f_Row(2,atom1) - dUdy
765 +       f_Row(3,atom1) = f_Row(3,atom1) - dUdz
766 +    
767 +       f_Col(1,atom2) = f_Col(1,atom2) + dUdx
768 +       f_Col(2,atom2) = f_Col(2,atom2) + dUdy
769 +       f_Col(3,atom2) = f_Col(3,atom2) + dUdz
770 +      
771 +       t_Row(1,atom2) = t_Row(1,atom1) + ul(2)*dUdu1z - ul(3)*dUdu1y
772 +       t_Row(2,atom2) = t_Row(2,atom1) + ul(3)*dUdu1x - ul(1)*dUdu1z
773 +       t_Row(3,atom2) = t_Row(3,atom1) + ul(1)*dUdu1y - ul(2)*dUdu1x
774 +  
775 + #else
776 +      
777 +       !!kdaily changed flx(gbatom) to f(1,atom1)
778 +       f(1,atom1) = f(1,atom1) + dUdx
779 +       f(2,atom1) = f(2,atom1) + dUdy
780 +       f(3,atom1) = f(3,atom1) + dUdz
781 +      
782 +       !!kdaily changed flx(ljatom) to f(2,atom2)
783 +       f(1,atom2) = f(1,atom2) - dUdx
784 +       f(2,atom2) = f(2,atom2) - dUdy
785 +       f(3,atom2) = f(3,atom2) - dUdz
786 +      
787 +       ! torques are cross products:
788 +       !!kdaily tlx(gbatom) to t(1, atom1)and ulx(gbatom) to u11(atom1)need mpi
789 +       t(1,atom2) = t(1,atom2) + ul(2)*dUduz - ul(3)*dUduy
790 +       t(2,atom2) = t(2,atom2) + ul(3)*dUdux - ul(1)*dUduz
791 +       t(3,atom2) = t(3,atom2) + ul(1)*dUduy - ul(2)*dUdux
792 +      
793 + #endif
794 +      
795 +       if (do_pot) then
796 + #ifdef IS_MPI
797 +          pot_row(atom1) = pot_row(atom1) + 2.0d0*eps*Rl26*sw
798 +          pot_col(atom2) = pot_col(atom2) + 2.0d0*eps*Rl26*sw
799 + #else
800 +          pot = pot + prefactor*eabf*Rl26*sw
801 + #endif
802 +       endif
803 +      
804 +       vpair = vpair + 4.0*epmu*Rl26
805 + #ifdef IS_MPI
806 +       id1 = AtomRowToGlobal(atom1)
807 +       id2 = AtomColToGlobal(atom2)
808 + #else
809 +       id1 = atom1
810 +       id2 = atom2
811 + #endif
812 +      
813 +       If (Molmembershiplist(Id1) .Ne. Molmembershiplist(Id2)) Then
814 +          
815 +          Fpair(1) = Fpair(1) + Dudx
816 +          Fpair(2) = Fpair(2) + Dudy
817 +          Fpair(3) = Fpair(3) + Dudz
818 +          
819 +       Endif
820 +      
821 +    else
822 +       !!need to do this all over but switch the gb and lj
823 +    endif
824      return
425 –  end subroutine do_gb_pair
825  
826 +  end subroutine do_gb_lj_pair
827 +
828 +  subroutine complete_GayBerne_FF(status)
829 +    integer :: nLJTYPES, nGayBerneTypes, ntypes, current, status, natypes
830 +    integer, pointer :: MatchList(:) => null ()
831 +    integer :: i
832 +    integer :: myATID
833 +    logical :: thisProperty
834 +    
835 +    if(.not.associated(LJMap%Ljtype)) then
836 +      
837 +       natypes = getSize(atypes)
838 +      
839 +       if(nAtypes == 0) then
840 +          status = -1
841 +          return
842 +       end if
843 +      
844 +       call getMatchingElementList(atypes, "is_LennardJones", .true., &
845 +            nLJTypes, MatchList)
846 +      
847 +       LJMap%nLJtypes = nLJTypes
848 +      
849 +       if(nLJTypes ==0) then
850 +          write(*,*)" you broke this thing kyle"
851 +          return
852 +       endif
853 +       allocate(LJMap%LJtype(nLJTypes))
854 +      
855 +       ntypes = getSize(atypes)
856 +      
857 +       allocate(LJMap%atidToLJtype(ntypes))
858 +    end if
859 +    
860 +    do i =1, ntypes
861 +      
862 +       myATID = getFirstMatchingElement(atypes, 'c_ident', i)
863 +       call getElementProperty(atypes, myATID, "is_LennardJones",thisProperty)
864 +      
865 +       if(thisProperty) then
866 +          current =  LJMap%currentLJtype+1      
867 +          LJMap%currentLJtype =  current          
868 +          
869 +          LJMap%atidToLJtype(myATID)    = current
870 +          LJMap%LJtype(current)%atid    = myATid
871 +          
872 +          LJMap%LJtype(current)%sigma   = getSigma(myATID)
873 +          LJMap%LJtype(current)%epsilon = getEpsilon(myATID)
874 +       endif
875 +      
876 +    enddo
877 +    gb_pair_initialized = .true.
878 +    
879 +  end subroutine complete_GayBerne_FF
880 +
881 +  subroutine destroyGayBerneTypes()
882 +
883 +    LJMap%Nljtypes =0
884 +    LJMap%currentLJtype=0
885 +    if(associated(LJMap%LJtype))then
886 +       deallocate(LJMap%LJtype)
887 +       LJMap%LJtype => null()
888 +    end if
889 +      
890 +    if(associated(LJMap%atidToLJType))then
891 +       deallocate(LJMap%atidToLJType)
892 +       LJMap%atidToLJType => null()
893 +    end if
894 +
895 + !!       deallocate(gayBerneLJMap)
896 +  
897 +    haveGayBerneLJMap = .false.
898 +  end subroutine destroyGayBerneTypes
899 +
900 +
901 +
902   end module gb_pair
903 +    

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines