OOPSE/libmdtools/calc_charge_charge.F90

module charge_charge
  
  use force_globals
  use definitions
  use atype_module
  use vector_class
  use simulation
  use status
#ifdef IS_MPI
  use mpiSimulation
#endif
  implicit none

  PRIVATE
  real(kind=dp), save :: ecr = 0.0_DP
  real(kind=dp), save :: rt  = 0.0_DP
  real(kind=dp), save :: dielect  = 0.0_DP
  real(kind=dp), save :: pre = 0.0_DP
  
  logical, save :: haveCutoffs = .false.
  logical, save :: haveChargeMap = .false.
  logical, save :: haveDielectric = .false.

  public::setCutoffsCharge
  public::do_charge_pair
  public::initialize_charge

  type :: ChargeList
     real(kind=DP) :: charge = 0.0_DP
  end type ChargeList

  type(ChargeList), dimension(:), allocatable :: ChargeMap

contains

  subroutine initialize_charge(this_dielect)
    real(kind=dp), intent(in) :: this_dielect
    
    dielect = this_dielect
    haveDielectric = .true.
    
    ! because setCutoffsCharge is called before initialize_charge
    ! we need to call it agin to make sure all of the precalculation
    ! value is correct. for the time being, just a quick hack
    call setCutoffsCharge(ecr, rt)
    return
  end subroutine initialize_charge
  
    
  subroutine setCutoffsCharge(this_ecr, this_rt)
    real(kind=dp), intent(in) :: this_ecr, this_rt
    ecr = this_ecr
    rt = this_rt

    ! pre converts from fundamental charge to kcal/mol
    pre = 332.06508_DP
    
    !if (.not.haveDielectric)then
    !   write(default_error,*) 'Dielect constant in charge module is not set!'
    !endif
    
    haveCutoffs = .true.
    
    return
  end subroutine setCutoffsCharge

  subroutine createChargeMap(status)
    integer :: nAtypes
    integer :: status
    integer :: i
    real (kind=DP) :: thisCharge
    logical :: thisProperty
    
    status = 0
    
    nAtypes = getSize(atypes)
    
    if (nAtypes == 0) then
       status = -1
       return
    end if
    
    if (.not. allocated(ChargeMap)) then
       allocate(ChargeMap(nAtypes))
    endif

    do i = 1, nAtypes

       call getElementProperty(atypes, i, "is_Charge", thisProperty)

       if (thisProperty) then
          call getElementProperty(atypes, i, "charge", thisCharge)
          ChargeMap(i)%charge = thisCharge
       endif
       
    end do
    
    haveChargeMap = .true.
    
  end subroutine createChargeMap
  
  
  subroutine do_charge_pair(atom1, atom2, d, rij, r2, dc, rcij, rc2, mfact, &
       pot, f, do_pot, do_stress, molecular_cutoffs)
    
    logical :: do_pot, do_stress, molecular_cutoffs

    integer atom1, atom2, me1, me2, id1, id2
    integer :: localError
    real(kind=dp) :: rij, r2, q1, q2, rcij, rc2
    real(kind=dp) :: drdx, drdy, drdz, dudr, fx, fy, fz
    real(kind=dp) :: taper, dtdr, vterm

    real( kind = dp ) :: pot
    real( kind = dp ), dimension(3) :: d, dc
    real( kind = dp ), dimension(3,nLocal) :: f
    real( kind = dp ), dimension(nLocal) :: mfact
    real( kind = dp ) :: theR, term1, term2
    real( kind = dp ) :: fxab, fyab, fzab, fxba, fyba, fzba
    
    if (.not. haveCutoffs) then
       write(default_error,*) 'charge-charge does not have cutoffs set!'
       return
    endif

    if (.not.haveChargeMap) then
       localError = 0
       call createChargeMap(localError)
       if ( localError .ne. 0 ) then
          call handleError("charge-charge", "ChargeMap creation failed!")
          return
       end if
    endif

#ifdef IS_MPI
    me1 = atid_Row(atom1)
    me2 = atid_Col(atom2)
#else
    me1 = atid(atom1)
    me2 = atid(atom2)
#endif

    q1 = ChargeMap(me1)%charge
    q2 = ChargeMap(me2)%charge

    if (molecular_cutoffs) then
       theR = rcij
    else
       theR = rij
    endif

    if (theR.le.ecr) then
       
       if (theR.lt.rt) then
          taper = 1.0d0
          dtdr = 0.0d0
       else
          taper = (ecr + 2.0d0*theR - 3.0d0*rt)*(ecr-theR)**2/ ((ecr-rt)**3)
          dtdr = 6.0d0*(theR*theR - theR*rt - theR*ecr +ecr*rt)/((ecr-rt)**3)
       endif

       vterm = pre * q1 * q2 / rij 

       if (molecular_cutoffs) then

          term1 = -vterm * (taper / rij)

          fx = term1 * d(1) / rij
          fy = term1 * d(2) / rij
          fz = term1 * d(3) / rij

          !term2 = vterm * dtdr * mfact(atom1)
          term2 = vterm * dtdr 

          fxab = term2 * dc(1) / rcij
          fyab = term2 * dc(2) / rcij
          fzab = term2 * dc(3) / rcij

          !term2 = vterm * dtdr * mfact(atom2)
          term2 = vterm * dtdr

          fxba = term2 * dc(1) / rcij
          fyba = term2 * dc(2) / rcij
          fzba = term2 * dc(3) / rcij

       else

          dudr =  vterm * (dtdr  - taper / rij)
          drdx = d(1) / rij
          drdy = d(2) / rij
          drdz = d(3) / rij
          
          fx = dudr * drdx
          fy = dudr * drdy
          fz = dudr * drdz          

       endif


#ifdef IS_MPI
       if (do_pot) then
          pot_Row(atom1) = pot_Row(atom1) + vterm*taper*0.5
          pot_Col(atom2) = pot_Col(atom2) + vterm*taper*0.5
       endif

       f_Row(1,atom1) = f_Row(1,atom1) + fx
       f_Row(2,atom1) = f_Row(2,atom1) + fy
       f_Row(3,atom1) = f_Row(3,atom1) + fz

       f_Col(1,atom2) = f_Col(1,atom2) - fx
       f_Col(2,atom2) = f_Col(2,atom2) - fy
       f_Col(3,atom2) = f_Col(3,atom2) - fz

       if (molecular_cutoffs) then
          f_Row(1,atom1) = f_Row(1,atom1) + fxab
          f_Row(2,atom1) = f_Row(2,atom1) + fyab
          f_Row(3,atom1) = f_Row(3,atom1) + fzab
          
          f_Col(1,atom2) = f_Col(1,atom2) - fxba
          f_Col(2,atom2) = f_Col(2,atom2) - fyba
          f_Col(3,atom2) = f_Col(3,atom2) - fzba
       endif

#else

       if (do_pot) pot = pot + vterm*taper

       f(1,atom1) = f(1,atom1) + fx
       f(2,atom1) = f(2,atom1) + fy
       f(3,atom1) = f(3,atom1) + fz

       f(1,atom2) = f(1,atom2) - fx
       f(2,atom2) = f(2,atom2) - fy
       f(3,atom2) = f(3,atom2) - fz

       if (molecular_cutoffs) then

          f(1,atom1) = f(1,atom1) + fxab
          f(2,atom1) = f(2,atom1) + fyab
          f(3,atom1) = f(3,atom1) + fzab
          
          f(1,atom2) = f(1,atom2) - fxba
          f(2,atom2) = f(2,atom2) - fyba
          f(3,atom2) = f(3,atom2) - fzba
       endif

#endif

       if (do_stress) then

#ifdef IS_MPI
          id1 = tagRow(atom1)
          id2 = tagColumn(atom2)
#else
          id1 = atom1
          id2 = atom2
#endif
 

          if (molMembershipList(id1) .ne. molMembershipList(id2)) then

             ! because the d vector is the rj - ri vector, and
             ! because fx, fy, fz are the force on atom i, we need a
             ! negative sign here:
             
             tau_Temp(1) = tau_Temp(1) - d(1) * fx
             tau_Temp(2) = tau_Temp(2) - d(1) * fy
             tau_Temp(3) = tau_Temp(3) - d(1) * fz
             tau_Temp(4) = tau_Temp(4) - d(2) * fx
             tau_Temp(5) = tau_Temp(5) - d(2) * fy
             tau_Temp(6) = tau_Temp(6) - d(2) * fz
             tau_Temp(8) = tau_Temp(8) - d(3) * fy
             tau_Temp(9) = tau_Temp(9) - d(3) * fz

             if (molecular_cutoffs) then
                
                tau_Temp(1) = tau_Temp(1) - d(1) * fxab
                tau_Temp(2) = tau_Temp(2) - d(1) * fyab
                tau_Temp(3) = tau_Temp(3) - d(1) * fzab
                tau_Temp(4) = tau_Temp(4) - d(2) * fxab
                tau_Temp(5) = tau_Temp(5) - d(2) * fyab
                tau_Temp(6) = tau_Temp(6) - d(2) * fzab
                tau_Temp(8) = tau_Temp(8) - d(3) * fyab
                tau_Temp(9) = tau_Temp(9) - d(3) * fzab
             endif
             
             virial_Temp = virial_Temp + &
                  (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
             
          endif
       endif
       
    endif
    return
  end subroutine do_charge_pair
  
end module charge_charge
Revision:	1144
Committed:	Sat May 1 18:52:38 2004 UTC (20 years, 2 months ago) by tim
File size:	7789 byte(s)
Log Message:	C++ pass groupList to fortran
#	User	Rev	Content
1	gezelter	940	module charge_charge
2
3			use force_globals
4			use definitions
5			use atype_module
6			use vector_class
7			use simulation
8			use status
9			#ifdef IS_MPI
10			use mpiSimulation
11			#endif
12			implicit none
13
14			PRIVATE
15	gezelter	945	real(kind=dp), save :: ecr = 0.0_DP
16	gezelter	940	real(kind=dp), save :: rt = 0.0_DP
17	tim	1132	real(kind=dp), save :: dielect = 0.0_DP
18	gezelter	940	real(kind=dp), save :: pre = 0.0_DP
19	tim	1132
20	gezelter	940	logical, save :: haveCutoffs = .false.
21			logical, save :: haveChargeMap = .false.
22	tim	1132	logical, save :: haveDielectric = .false.
23	gezelter	940
24	gezelter	945	public::setCutoffsCharge
25	gezelter	940	public::do_charge_pair
26	tim	1132	public::initialize_charge
27	gezelter	940
28			type :: ChargeList
29			real(kind=DP) :: charge = 0.0_DP
30			end type ChargeList
31
32			type(ChargeList), dimension(:), allocatable :: ChargeMap
33
34			contains
35	tim	1132
36			subroutine initialize_charge(this_dielect)
37			real(kind=dp), intent(in) :: this_dielect
38	gezelter	940
39	tim	1132	dielect = this_dielect
40			haveDielectric = .true.
41
42			! because setCutoffsCharge is called before initialize_charge
43			! we need to call it agin to make sure all of the precalculation
44			! value is correct. for the time being, just a quick hack
45			call setCutoffsCharge(ecr, rt)
46			return
47			end subroutine initialize_charge
48
49
50	gezelter	945	subroutine setCutoffsCharge(this_ecr, this_rt)
51			real(kind=dp), intent(in) :: this_ecr, this_rt
52			ecr = this_ecr
53	gezelter	940	rt = this_rt
54
55			! pre converts from fundamental charge to kcal/mol
56			pre = 332.06508_DP
57	tim	1132
58			!if (.not.haveDielectric)then
59			! write(default_error,*) 'Dielect constant in charge module is not set!'
60			!endif
61
62	gezelter	940	haveCutoffs = .true.
63
64			return
65	gezelter	945	end subroutine setCutoffsCharge
66	gezelter	940
67			subroutine createChargeMap(status)
68			integer :: nAtypes
69			integer :: status
70			integer :: i
71			real (kind=DP) :: thisCharge
72			logical :: thisProperty
73
74			status = 0
75
76			nAtypes = getSize(atypes)
77
78			if (nAtypes == 0) then
79			status = -1
80			return
81			end if
82
83			if (.not. allocated(ChargeMap)) then
84			allocate(ChargeMap(nAtypes))
85			endif
86
87			do i = 1, nAtypes
88
89			call getElementProperty(atypes, i, "is_Charge", thisProperty)
90
91			if (thisProperty) then
92			call getElementProperty(atypes, i, "charge", thisCharge)
93			ChargeMap(i)%charge = thisCharge
94			endif
95
96			end do
97
98			haveChargeMap = .true.
99
100			end subroutine createChargeMap
101
102
103	gezelter	1138	subroutine do_charge_pair(atom1, atom2, d, rij, r2, dc, rcij, rc2, mfact, &
104			pot, f, do_pot, do_stress, molecular_cutoffs)
105	gezelter	940
106	gezelter	1138	logical :: do_pot, do_stress, molecular_cutoffs
107	gezelter	940
108			integer atom1, atom2, me1, me2, id1, id2
109			integer :: localError
110	gezelter	1138	real(kind=dp) :: rij, r2, q1, q2, rcij, rc2
111	gezelter	940	real(kind=dp) :: drdx, drdy, drdz, dudr, fx, fy, fz
112			real(kind=dp) :: taper, dtdr, vterm
113
114			real( kind = dp ) :: pot
115	gezelter	1138	real( kind = dp ), dimension(3) :: d, dc
116	gezelter	940	real( kind = dp ), dimension(3,nLocal) :: f
117	gezelter	1138	real( kind = dp ), dimension(nLocal) :: mfact
118	tim	1143	real( kind = dp ) :: theR, term1, term2
119	gezelter	1138	real( kind = dp ) :: fxab, fyab, fzab, fxba, fyba, fzba
120	gezelter	940
121			if (.not. haveCutoffs) then
122			write(default_error,*) 'charge-charge does not have cutoffs set!'
123			return
124			endif
125
126			if (.not.haveChargeMap) then
127			localError = 0
128			call createChargeMap(localError)
129			if ( localError .ne. 0 ) then
130			call handleError("charge-charge", "ChargeMap creation failed!")
131			return
132			end if
133			endif
134
135			#ifdef IS_MPI
136			me1 = atid_Row(atom1)
137			me2 = atid_Col(atom2)
138			#else
139			me1 = atid(atom1)
140			me2 = atid(atom2)
141			#endif
142
143			q1 = ChargeMap(me1)%charge
144			q2 = ChargeMap(me2)%charge
145	gezelter	1138
146			if (molecular_cutoffs) then
147			theR = rcij
148			else
149			theR = rij
150			endif
151
152			if (theR.le.ecr) then
153	gezelter	940
154	gezelter	1138	if (theR.lt.rt) then
155	gezelter	940	taper = 1.0d0
156			dtdr = 0.0d0
157			else
158	gezelter	1138	taper = (ecr + 2.0d0theR - 3.0d0rt)(ecr-theR)2/ ((ecr-rt)*3)
159			dtdr = 6.0d0(theRtheR - theRrt - theRecr +ecrrt)/((ecr-rt)*3)
160	gezelter	940	endif
161
162	gezelter	1138	vterm = pre * q1 * q2 / rij
163	gezelter	1133
164	gezelter	1138	if (molecular_cutoffs) then
165	gezelter	940
166	tim	1143	term1 = -vterm * (taper / rij)
167	gezelter	940
168	tim	1143	fx = term1 * d(1) / rij
169			fy = term1 * d(2) / rij
170			fz = term1 * d(3) / rij
171
172	tim	1144	!term2 = vterm * dtdr * mfact(atom1)
173			term2 = vterm * dtdr
174	tim	1113
175	gezelter	1138	fxab = term2 * dc(1) / rcij
176			fyab = term2 * dc(2) / rcij
177			fzab = term2 * dc(3) / rcij
178
179	tim	1144	!term2 = vterm * dtdr * mfact(atom2)
180			term2 = vterm * dtdr
181	gezelter	1138
182			fxba = term2 * dc(1) / rcij
183			fyba = term2 * dc(2) / rcij
184			fzba = term2 * dc(3) / rcij
185
186			else
187
188			dudr = vterm * (dtdr - taper / rij)
189			drdx = d(1) / rij
190			drdy = d(2) / rij
191			drdz = d(3) / rij
192
193			fx = dudr * drdx
194			fy = dudr * drdy
195			fz = dudr * drdz
196
197			endif
198
199
200	gezelter	940	#ifdef IS_MPI
201			if (do_pot) then
202			pot_Row(atom1) = pot_Row(atom1) + vtermtaper0.5
203			pot_Col(atom2) = pot_Col(atom2) + vtermtaper0.5
204			endif
205
206			f_Row(1,atom1) = f_Row(1,atom1) + fx
207			f_Row(2,atom1) = f_Row(2,atom1) + fy
208			f_Row(3,atom1) = f_Row(3,atom1) + fz
209
210			f_Col(1,atom2) = f_Col(1,atom2) - fx
211			f_Col(2,atom2) = f_Col(2,atom2) - fy
212			f_Col(3,atom2) = f_Col(3,atom2) - fz
213
214	gezelter	1138	if (molecular_cutoffs) then
215			f_Row(1,atom1) = f_Row(1,atom1) + fxab
216			f_Row(2,atom1) = f_Row(2,atom1) + fyab
217			f_Row(3,atom1) = f_Row(3,atom1) + fzab
218
219			f_Col(1,atom2) = f_Col(1,atom2) - fxba
220			f_Col(2,atom2) = f_Col(2,atom2) - fyba
221			f_Col(3,atom2) = f_Col(3,atom2) - fzba
222			endif
223
224	gezelter	940	#else
225
226			if (do_pot) pot = pot + vterm*taper
227
228			f(1,atom1) = f(1,atom1) + fx
229			f(2,atom1) = f(2,atom1) + fy
230			f(3,atom1) = f(3,atom1) + fz
231
232			f(1,atom2) = f(1,atom2) - fx
233			f(2,atom2) = f(2,atom2) - fy
234			f(3,atom2) = f(3,atom2) - fz
235
236	gezelter	1138	if (molecular_cutoffs) then
237	tim	1144
238	gezelter	1138	f(1,atom1) = f(1,atom1) + fxab
239			f(2,atom1) = f(2,atom1) + fyab
240			f(3,atom1) = f(3,atom1) + fzab
241
242			f(1,atom2) = f(1,atom2) - fxba
243			f(2,atom2) = f(2,atom2) - fyba
244			f(3,atom2) = f(3,atom2) - fzba
245			endif
246
247	gezelter	940	#endif
248
249			if (do_stress) then
250
251			#ifdef IS_MPI
252			id1 = tagRow(atom1)
253			id2 = tagColumn(atom2)
254			#else
255			id1 = atom1
256			id2 = atom2
257			#endif
258	tim	1131
259	gezelter	940
260			if (molMembershipList(id1) .ne. molMembershipList(id2)) then
261	tim	1131
262	gezelter	940	! because the d vector is the rj - ri vector, and
263			! because fx, fy, fz are the force on atom i, we need a
264			! negative sign here:
265
266			tau_Temp(1) = tau_Temp(1) - d(1) * fx
267			tau_Temp(2) = tau_Temp(2) - d(1) * fy
268			tau_Temp(3) = tau_Temp(3) - d(1) * fz
269			tau_Temp(4) = tau_Temp(4) - d(2) * fx
270			tau_Temp(5) = tau_Temp(5) - d(2) * fy
271			tau_Temp(6) = tau_Temp(6) - d(2) * fz
272			tau_Temp(8) = tau_Temp(8) - d(3) * fy
273			tau_Temp(9) = tau_Temp(9) - d(3) * fz
274	gezelter	1138
275			if (molecular_cutoffs) then
276
277			tau_Temp(1) = tau_Temp(1) - d(1) * fxab
278			tau_Temp(2) = tau_Temp(2) - d(1) * fyab
279			tau_Temp(3) = tau_Temp(3) - d(1) * fzab
280			tau_Temp(4) = tau_Temp(4) - d(2) * fxab
281			tau_Temp(5) = tau_Temp(5) - d(2) * fyab
282			tau_Temp(6) = tau_Temp(6) - d(2) * fzab
283			tau_Temp(8) = tau_Temp(8) - d(3) * fyab
284			tau_Temp(9) = tau_Temp(9) - d(3) * fzab
285			endif
286	gezelter	940
287			virial_Temp = virial_Temp + &
288			(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
289
290			endif
291			endif
292
293			endif
294			return
295			end subroutine do_charge_pair
296
297			end module charge_charge