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

Comparing trunk/OOPSE-3.0/src/UseTheForce/DarkSide/electrostatic.F90 (file contents):
Revision 2189 by chuckv, Wed Apr 13 20:36:45 2005 UTC vs.
Revision 2204 by gezelter, Fri Apr 15 22:04:00 2005 UTC

# Line 40 | Line 40 | module electrostatic_module
40   !!
41  
42   module electrostatic_module
43 <  
43 >
44    use force_globals
45    use definitions
46    use atype_module
# Line 97 | Line 97 | contains
97  
98    subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, &
99         is_SplitDipole, is_Quadrupole, status)
100 <    
100 >
101      integer, intent(in) :: c_ident
102      logical, intent(in) :: is_Charge
103      logical, intent(in) :: is_Dipole
# Line 108 | Line 108 | contains
108  
109      status = 0
110      myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
111 <    
111 >
112      !! Be simple-minded and assume that we need an ElectrostaticMap that
113      !! is the same size as the total number of atom types
114  
115      if (.not.allocated(ElectrostaticMap)) then
116 <      
116 >
117         nAtypes = getSize(atypes)
118 <    
118 >
119         if (nAtypes == 0) then
120            status = -1
121            return
122         end if
123 <      
123 >
124         if (.not. allocated(ElectrostaticMap)) then
125            allocate(ElectrostaticMap(nAtypes))
126         endif
127 <      
127 >
128      end if
129  
130      if (myATID .gt. size(ElectrostaticMap)) then
131         status = -1
132         return
133      endif
134 <    
134 >
135      ! set the values for ElectrostaticMap for this atom type:
136  
137      ElectrostaticMap(myATID)%c_ident = c_ident
# Line 139 | Line 139 | contains
139      ElectrostaticMap(myATID)%is_Dipole = is_Dipole
140      ElectrostaticMap(myATID)%is_SplitDipole = is_SplitDipole
141      ElectrostaticMap(myATID)%is_Quadrupole = is_Quadrupole
142 <    
142 >
143    end subroutine newElectrostaticType
144  
145    subroutine setCharge(c_ident, charge, status)
# Line 167 | Line 167 | contains
167         call handleError("electrostatic", "Attempt to setCharge of an atom type that is not a charge!")
168         status = -1
169         return
170 <    endif      
170 >    endif
171  
172      ElectrostaticMap(myATID)%charge = charge
173    end subroutine setCharge
# Line 258 | Line 258 | contains
258         status = -1
259         return
260      endif
261 <    
261 >
262      do i = 1, 3
263 <          ElectrostaticMap(myATID)%quadrupole_moments(i) = &
264 <               quadrupole_moments(i)
265 <       enddo
263 >       ElectrostaticMap(myATID)%quadrupole_moments(i) = &
264 >            quadrupole_moments(i)
265 >    enddo
266  
267    end subroutine setQuadrupoleMoments
268  
269 <  
269 >
270    function getCharge(atid) result (c)
271      integer, intent(in) :: atid
272      integer :: localError
273      real(kind=dp) :: c
274 <    
274 >
275      if (.not.allocated(ElectrostaticMap)) then
276         call handleError("electrostatic", "no ElectrostaticMap was present before first call of getCharge!")
277         return
278      end if
279 <    
279 >
280      if (.not.ElectrostaticMap(atid)%is_Charge) then
281         call handleError("electrostatic", "getCharge was called for an atom type that isn't a charge!")
282         return
283      endif
284 <    
284 >
285      c = ElectrostaticMap(atid)%charge
286    end function getCharge
287  
# Line 289 | Line 289 | contains
289      integer, intent(in) :: atid
290      integer :: localError
291      real(kind=dp) :: dm
292 <    
292 >
293      if (.not.allocated(ElectrostaticMap)) then
294         call handleError("electrostatic", "no ElectrostaticMap was present before first call of getDipoleMoment!")
295         return
296      end if
297 <    
297 >
298      if (.not.ElectrostaticMap(atid)%is_Dipole) then
299         call handleError("electrostatic", "getDipoleMoment was called for an atom type that isn't a dipole!")
300         return
301      endif
302 <    
302 >
303      dm = ElectrostaticMap(atid)%dipole_moment
304    end function getDipoleMoment
305  
306    subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
307         vpair, fpair, pot, eFrame, f, t, do_pot)
308 <    
308 >
309      logical, intent(in) :: do_pot
310 <    
310 >
311      integer, intent(in) :: atom1, atom2
312      integer :: localError
313  
# Line 320 | Line 320 | contains
320      real( kind = dp ), dimension(9,nLocal) :: eFrame
321      real( kind = dp ), dimension(3,nLocal) :: f
322      real( kind = dp ), dimension(3,nLocal) :: t
323 <    
323 >
324      real (kind = dp), dimension(3) :: ux_i, uy_i, uz_i
325      real (kind = dp), dimension(3) :: ux_j, uy_j, uz_j
326      real (kind = dp), dimension(3) :: dudux_i, duduy_i, duduz_i
# Line 378 | Line 378 | contains
378      if (i_is_Charge) then
379         q_i = ElectrostaticMap(me1)%charge      
380      endif
381 <    
381 >
382      if (i_is_Dipole) then
383         mu_i = ElectrostaticMap(me1)%dipole_moment
384   #ifdef IS_MPI
# Line 395 | Line 395 | contains
395         if (i_is_SplitDipole) then
396            d_i = ElectrostaticMap(me1)%split_dipole_distance
397         endif
398 <      
398 >
399      endif
400  
401      if (i_is_Quadrupole) then
# Line 432 | Line 432 | contains
432      if (j_is_Charge) then
433         q_j = ElectrostaticMap(me2)%charge      
434      endif
435 <    
435 >
436      if (j_is_Dipole) then
437         mu_j = ElectrostaticMap(me2)%dipole_moment
438   #ifdef IS_MPI
# Line 497 | Line 497 | contains
497      if (i_is_Charge) then
498  
499         if (j_is_Charge) then
500 <          
500 >
501            vterm = pre11 * q_i * q_j * riji
502            vpair = vpair + vterm
503            epot = epot + sw*vterm
# Line 507 | Line 507 | contains
507            dudx = dudx + dudr * xhat
508            dudy = dudy + dudr * yhat
509            dudz = dudz + dudr * zhat
510 <      
510 >
511         endif
512  
513         if (j_is_Dipole) then
# Line 524 | Line 524 | contains
524            ri2 = ri * ri
525            ri3 = ri2 * ri
526            sc2 = scale * scale
527 <            
527 >
528            pref = pre12 * q_i * mu_j
529            vterm = - pref * ct_j * ri2 * scale
530            vpair = vpair + vterm
# Line 541 | Line 541 | contains
541            duduz_j(1) = duduz_j(1) - pref * sw * ri2 * xhat * scale
542            duduz_j(2) = duduz_j(2) - pref * sw * ri2 * yhat * scale
543            duduz_j(3) = duduz_j(3) - pref * sw * ri2 * zhat * scale
544 <          
544 >
545         endif
546  
547         if (j_is_Quadrupole) then
# Line 572 | Line 572 | contains
572                 qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
573                 qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
574                 qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
575 <          
575 >
576            dudux_j(1) = dudux_j(1) + pref * sw * ri3 * (qxx_j*6.0_dp*cx_j*xhat)
577            dudux_j(2) = dudux_j(2) + pref * sw * ri3 * (qxx_j*6.0_dp*cx_j*yhat)
578            dudux_j(3) = dudux_j(3) + pref * sw * ri3 * (qxx_j*6.0_dp*cx_j*zhat)
# Line 587 | Line 587 | contains
587         endif
588  
589      endif
590 <  
590 >
591      if (i_is_Dipole) then
592 <      
592 >
593         if (j_is_Charge) then
594  
595            if (i_is_SplitDipole) then
# Line 604 | Line 604 | contains
604            ri2 = ri * ri
605            ri3 = ri2 * ri
606            sc2 = scale * scale
607 <            
607 >
608            pref = pre12 * q_j * mu_i
609            vterm = pref * ct_i * ri2 * scale
610            vpair = vpair + vterm
# Line 651 | Line 651 | contains
651            vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
652            vpair = vpair + vterm
653            epot = epot + sw * vterm
654 <          
654 >
655            a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
656  
657            dudx=dudx+pref*sw*3.0d0*ri4*scale*(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
# Line 671 | Line 671 | contains
671  
672      if (i_is_Quadrupole) then
673         if (j_is_Charge) then
674 <          
674 >
675            ri2 = riji * riji
676            ri3 = ri2 * riji
677            ri4 = ri2 * ri2
678            cx2 = cx_i * cx_i
679            cy2 = cy_i * cy_i
680            cz2 = cz_i * cz_i
681 <          
681 >
682            pref = pre14 * q_j / 3.0_dp
683            vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
684                 qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
685                 qzz_i * (3.0_dp*cz2 - 1.0_dp))
686            vpair = vpair + vterm
687            epot = epot + sw * vterm
688 <          
688 >
689            dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + pref * sw * ri4 * ( &
690                 qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
691                 qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
# Line 698 | Line 698 | contains
698                 qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
699                 qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
700                 qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
701 <          
701 >
702            dudux_i(1) = dudux_i(1) + pref * sw * ri3 * (qxx_i*6.0_dp*cx_i*xhat)
703            dudux_i(2) = dudux_i(2) + pref * sw * ri3 * (qxx_i*6.0_dp*cx_i*yhat)
704            dudux_i(3) = dudux_i(3) + pref * sw * ri3 * (qxx_i*6.0_dp*cx_i*zhat)
705 <          
705 >
706            duduy_i(1) = duduy_i(1) + pref * sw * ri3 * (qyy_i*6.0_dp*cy_i*xhat)
707            duduy_i(2) = duduy_i(2) + pref * sw * ri3 * (qyy_i*6.0_dp*cy_i*yhat)
708            duduy_i(3) = duduy_i(3) + pref * sw * ri3 * (qyy_i*6.0_dp*cy_i*zhat)
709 <          
709 >
710            duduz_i(1) = duduz_i(1) + pref * sw * ri3 * (qzz_i*6.0_dp*cz_i*xhat)
711            duduz_i(2) = duduz_i(2) + pref * sw * ri3 * (qzz_i*6.0_dp*cz_i*yhat)
712            duduz_i(3) = duduz_i(3) + pref * sw * ri3 * (qzz_i*6.0_dp*cz_i*zhat)
713         endif
714      endif
715 <      
716 <    
715 >
716 >
717      if (do_pot) then
718   #ifdef IS_MPI
719         pot_row(atom1) = pot_row(atom1) + 0.5d0*epot
# Line 722 | Line 722 | contains
722         pot = pot + epot
723   #endif
724      endif
725 <        
725 >
726   #ifdef IS_MPI
727      f_Row(1,atom1) = f_Row(1,atom1) + dudx
728      f_Row(2,atom1) = f_Row(2,atom1) + dudy
729      f_Row(3,atom1) = f_Row(3,atom1) + dudz
730 <    
730 >
731      f_Col(1,atom2) = f_Col(1,atom2) - dudx
732      f_Col(2,atom2) = f_Col(2,atom2) - dudy
733      f_Col(3,atom2) = f_Col(3,atom2) - dudz
734 <    
734 >
735      if (i_is_Dipole .or. i_is_Quadrupole) then
736         t_Row(1,atom1)=t_Row(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
737         t_Row(2,atom1)=t_Row(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
# Line 766 | Line 766 | contains
766      f(1,atom1) = f(1,atom1) + dudx
767      f(2,atom1) = f(2,atom1) + dudy
768      f(3,atom1) = f(3,atom1) + dudz
769 <    
769 >
770      f(1,atom2) = f(1,atom2) - dudx
771      f(2,atom2) = f(2,atom2) - dudy
772      f(3,atom2) = f(3,atom2) - dudz
773 <    
773 >
774      if (i_is_Dipole .or. i_is_Quadrupole) then
775         t(1,atom1)=t(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
776         t(2,atom1)=t(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
# Line 802 | Line 802 | contains
802      endif
803  
804   #endif
805 <    
805 >
806   #ifdef IS_MPI
807      id1 = AtomRowToGlobal(atom1)
808      id2 = AtomColToGlobal(atom2)
# Line 812 | Line 812 | contains
812   #endif
813  
814      if (molMembershipList(id1) .ne. molMembershipList(id2)) then
815 <      
815 >
816         fpair(1) = fpair(1) + dudx
817         fpair(2) = fpair(2) + dudy
818         fpair(3) = fpair(3) + dudz
# Line 821 | Line 821 | contains
821  
822      return
823    end subroutine doElectrostaticPair
824 –  
824  
825 +
826    subroutine destroyElectrostaticTypes()
827 –    
828 –   if(allocated(ElectrostaticMap)) deallocate(ElectrostaticMap)
827  
828 +    if(allocated(ElectrostaticMap)) deallocate(ElectrostaticMap)
829 +
830    end subroutine destroyElectrostaticTypes
831  
832   end module electrostatic_module

Diff Legend

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