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 1930 by gezelter, Wed Jan 12 22:41:40 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    
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      
# Line 157 | Line 333 | contains
333      ul2(1) = A(3,atom2)
334      ul2(2) = A(6,atom2)
335      ul2(3) = A(9,atom2)
336 +
337   #endif
338      
339      dru1dx = ul1(1)
# Line 165 | Line 342 | contains
342      dru2dy = ul2(2)
343      dru1dz = ul1(3)
344      dru2dz = ul2(3)
345 <    
345 >        
346 >
347 >
348      drdx = d(1) / r
349      drdy = d(2) / r
350      drdz = d(3) / r
# 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
447      dBigRdu2x = dgdu2x*gfact
448      dBigRdu2y = dgdu2y*gfact
449      dBigRdu2z = dgdu2z*gfact
450 <  
450 >
451      ! Now, we must do it again for g(ChiPrime) and dgpdx
452  
453      line1a = dotsum/opXpdot
# Line 286 | Line 466 | contains
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)
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 <    
475 >
476      term2a = -line3a/opXpdot
477      term2b =  line3b/omXpdot
478      
# Line 316 | Line 496 | contains
496      gmu = gp**gb_mu
497      gpi = 1.0d0 / gp
498      gmum = gmu*gpi
319 –  
320 –    ! write(*,*) atom1, atom2, Chi, u1dotu2
321 –    curlyE = 1.0d0/dsqrt(1.0d0 - Chi*Chi*u1dotu2*u1dotu2)
499  
500 +    curlyE = 1.0d0/dsqrt(1.0d0 - Chi*Chi*u1dotu2*u1dotu2)
501 + !!$
502 + !!$    dcE = -(curlyE**3)*Chi*Chi*u1dotu2
503      dcE = (curlyE**3)*Chi*Chi*u1dotu2
504 <  
504 >
505      dcEdu1x = dcE*ul2(1)
506      dcEdu1y = dcE*ul2(2)
507      dcEdu1z = dcE*ul2(3)
# Line 370 | Line 550 | contains
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
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(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
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
# Line 386 | Line 566 | contains
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
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(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
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 424 | Line 604 | end module gb_pair
604      return
605    end subroutine do_gb_pair
606  
607 < end module gb_pair
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 <  subroutine set_gb_pair_params(sigma, l2b_ratio, eps, eps_ratio, mu, nu)
661 <    use definitions, ONLY : dp
662 <    use gb_pair, ONLY : module_set_gb_pair_params => set_gb_pair_params
663 <    real( kind = dp ), intent(inout) :: sigma, l2b_ratio, eps, eps_ratio
664 <    real( kind = dp ), intent(inout) :: mu, nu
665 <    call module_set_gb_pair_params(sigma, l2b_ratio, eps, eps_ratio, mu, nu)
666 < end subroutine set_gb_pair_params
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
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