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, fxij, fyij, fzij
    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)
         
          fxij = term1 * d(1) / rij
          fyij = term1 * d(2) / rij
          fzij = term1 * d(3) / rij

          term2 = vterm * dtdr * mfact(atom1)

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

          term2 = vterm * dtdr * mfact(atom2)

          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:	1138
Committed:	Wed Apr 28 21:39:22 2004 UTC (20 years, 2 months ago) by gezelter
File size:	7755 byte(s)
Log Message:	work on molecular cutoffs
#	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			real( kind = dp ) :: theR, term1, term2, fxij, fyij, fzij
119			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	gezelter	1138	term1 = vterm * (taper / rij)
167
168			fxij = term1 * d(1) / rij
169			fyij = term1 * d(2) / rij
170			fzij = term1 * d(3) / rij
171	gezelter	940
172	gezelter	1138	term2 = vterm * dtdr * mfact(atom1)
173	tim	1113
174	gezelter	1138	fxab = term2 * dc(1) / rcij
175			fyab = term2 * dc(2) / rcij
176			fzab = term2 * dc(3) / rcij
177
178			term2 = vterm * dtdr * mfact(atom2)
179
180			fxba = term2 * dc(1) / rcij
181			fyba = term2 * dc(2) / rcij
182			fzba = term2 * dc(3) / rcij
183
184			else
185
186			dudr = vterm * (dtdr - taper / rij)
187			drdx = d(1) / rij
188			drdy = d(2) / rij
189			drdz = d(3) / rij
190
191			fx = dudr * drdx
192			fy = dudr * drdy
193			fz = dudr * drdz
194
195			endif
196
197
198	gezelter	940	#ifdef IS_MPI
199			if (do_pot) then
200			pot_Row(atom1) = pot_Row(atom1) + vtermtaper0.5
201			pot_Col(atom2) = pot_Col(atom2) + vtermtaper0.5
202			endif
203
204			f_Row(1,atom1) = f_Row(1,atom1) + fx
205			f_Row(2,atom1) = f_Row(2,atom1) + fy
206			f_Row(3,atom1) = f_Row(3,atom1) + fz
207
208			f_Col(1,atom2) = f_Col(1,atom2) - fx
209			f_Col(2,atom2) = f_Col(2,atom2) - fy
210			f_Col(3,atom2) = f_Col(3,atom2) - fz
211
212	gezelter	1138	if (molecular_cutoffs) then
213			f_Row(1,atom1) = f_Row(1,atom1) + fxab
214			f_Row(2,atom1) = f_Row(2,atom1) + fyab
215			f_Row(3,atom1) = f_Row(3,atom1) + fzab
216
217			f_Col(1,atom2) = f_Col(1,atom2) - fxba
218			f_Col(2,atom2) = f_Col(2,atom2) - fyba
219			f_Col(3,atom2) = f_Col(3,atom2) - fzba
220			endif
221
222	gezelter	940	#else
223
224			if (do_pot) pot = pot + vterm*taper
225
226			f(1,atom1) = f(1,atom1) + fx
227			f(2,atom1) = f(2,atom1) + fy
228			f(3,atom1) = f(3,atom1) + fz
229
230			f(1,atom2) = f(1,atom2) - fx
231			f(2,atom2) = f(2,atom2) - fy
232			f(3,atom2) = f(3,atom2) - fz
233
234	gezelter	1138	if (molecular_cutoffs) then
235			f(1,atom1) = f(1,atom1) + fxab
236			f(2,atom1) = f(2,atom1) + fyab
237			f(3,atom1) = f(3,atom1) + fzab
238
239			f(1,atom2) = f(1,atom2) - fxba
240			f(2,atom2) = f(2,atom2) - fyba
241			f(3,atom2) = f(3,atom2) - fzba
242			endif
243
244	gezelter	940	#endif
245
246			if (do_stress) then
247
248			#ifdef IS_MPI
249			id1 = tagRow(atom1)
250			id2 = tagColumn(atom2)
251			#else
252			id1 = atom1
253			id2 = atom2
254			#endif
255	tim	1131
256	gezelter	940
257			if (molMembershipList(id1) .ne. molMembershipList(id2)) then
258	tim	1131
259	gezelter	940	! because the d vector is the rj - ri vector, and
260			! because fx, fy, fz are the force on atom i, we need a
261			! negative sign here:
262
263			tau_Temp(1) = tau_Temp(1) - d(1) * fx
264			tau_Temp(2) = tau_Temp(2) - d(1) * fy
265			tau_Temp(3) = tau_Temp(3) - d(1) * fz
266			tau_Temp(4) = tau_Temp(4) - d(2) * fx
267			tau_Temp(5) = tau_Temp(5) - d(2) * fy
268			tau_Temp(6) = tau_Temp(6) - d(2) * fz
269			tau_Temp(8) = tau_Temp(8) - d(3) * fy
270			tau_Temp(9) = tau_Temp(9) - d(3) * fz
271	gezelter	1138
272			if (molecular_cutoffs) then
273
274			tau_Temp(1) = tau_Temp(1) - d(1) * fxab
275			tau_Temp(2) = tau_Temp(2) - d(1) * fyab
276			tau_Temp(3) = tau_Temp(3) - d(1) * fzab
277			tau_Temp(4) = tau_Temp(4) - d(2) * fxab
278			tau_Temp(5) = tau_Temp(5) - d(2) * fyab
279			tau_Temp(6) = tau_Temp(6) - d(2) * fzab
280			tau_Temp(8) = tau_Temp(8) - d(3) * fyab
281			tau_Temp(9) = tau_Temp(9) - d(3) * fzab
282			endif
283	gezelter	940
284			virial_Temp = virial_Temp + &
285			(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
286
287			endif
288			endif
289
290			endif
291			return
292			end subroutine do_charge_pair
293
294			end module charge_charge