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 2361 by gezelter, Wed Oct 12 21:00:50 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
54 < #define __FORTRAN90
55 < #include "UseTheForce/DarkSide/fInteractionMap.h"
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 76 | 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)
# Line 168 | 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 176 | 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 231 | Line 399 | contains
399      dgdy = term1y + line3y
400      dgdz = term1z + line3z
401  
402 <    term1u1x = (line1a+line2a)*dru1dx
403 <    term1u1y = (line1a+line2a)*dru1dy
404 <    term1u1z = (line1a+line2a)*dru1dz
405 <    term1u2x = (line1a-line2a)*dru2dx
406 <    term1u2y = (line1a-line2a)*dru2dy
407 <    term1u2z = (line1a-line2a)*dru2dz
402 >    term1u1x = 2.0d0*(line1a+line2a)*d(1)
403 >    term1u1y = 2.0d0*(line1a+line2a)*d(2)
404 >    term1u1z = 2.0d0*(line1a+line2a)*d(3)
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      
409      term2a = -line3a/opXdot
410      term2b =  line3b/omXdot
# Line 298 | Line 466 | contains
466      dgpdy = term1y + line3y
467      dgpdz = term1z + line3z
468      
469 <    term1u1x = (line1a+line2a)*dru1dx
470 <    term1u1y = (line1a+line2a)*dru1dy
471 <    term1u1z = (line1a+line2a)*dru1dz
472 <    term1u2x = (line1a-line2a)*dru2dx
473 <    term1u2y = (line1a-line2a)*dru2dy
474 <    term1u2z = (line1a-line2a)*dru2dz
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
# Line 330 | Line 498 | contains
498      gmum = gmu*gpi
499  
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 <  
504 >
505      dcEdu1x = dcE*ul2(1)
506      dcEdu1y = dcE*ul2(2)
507      dcEdu1z = dcE*ul2(3)
# Line 381 | 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 397 | 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(VDW_POT,atom1) = pot_row(VDW_POT,atom1) + 2.0d0*eps*R126*sw
581 <       pot_col(VDW_POT,atom2) = pot_col(VDW_POT,atom2) + 2.0d0*eps*R126*sw
580 >       pot_row(atom1) = pot_row(atom1) + 2.0d0*eps*R126*sw
581 >       pot_col(atom2) = pot_col(atom2) + 2.0d0*eps*R126*sw
582   #else
583         pot = pot + 4.0*eps*R126*sw
584   #endif
# Line 435 | Line 604 | end module gb_pair
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
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