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, sw, vpair, &
       pot, f, do_pot, do_stress)
    
    logical :: do_pot, do_stress

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

    real( kind = dp ) :: pot
    real( kind = dp ), dimension(3) :: d
    real( kind = dp ), dimension(3,nLocal) :: f

    
    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

    vterm = sw * pre * q1 * q2 / rij 

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

    vpair = vpair + vterm
    
#ifdef IS_MPI
    if (do_pot) then
       pot_Row(atom1) = pot_Row(atom1) + vterm*0.5
       pot_Col(atom2) = pot_Col(atom2) + vterm*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
    
#else

    if (do_pot) pot = pot + vterm
    
    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
        
#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
          
          
          virial_Temp = virial_Temp + &
               (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
          
       endif
    endif
    
    
    return
  end subroutine do_charge_pair
  
end module charge_charge
Revision:	1150
Committed:	Fri May 7 21:35:05 2004 UTC (20 years, 2 months ago) by gezelter
File size:	5299 byte(s)
Log Message:	Many changes to get group-based cutoffs to work
#	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	1150	subroutine do_charge_pair(atom1, atom2, d, rij, r2, sw, vpair, &
104			pot, f, do_pot, do_stress)
105	gezelter	940
106	gezelter	1150	logical :: do_pot, do_stress
107	gezelter	940
108			integer atom1, atom2, me1, me2, id1, id2
109			integer :: localError
110	gezelter	1150	real(kind=dp) :: rij, r2, q1, q2, sw, vpair
111	gezelter	940	real(kind=dp) :: drdx, drdy, drdz, dudr, fx, fy, fz
112	gezelter	1150	real(kind=dp) :: vterm
113	gezelter	940
114			real( kind = dp ) :: pot
115	gezelter	1150	real( kind = dp ), dimension(3) :: d
116	gezelter	940	real( kind = dp ), dimension(3,nLocal) :: f
117	gezelter	1150
118	gezelter	940
119			if (.not. haveCutoffs) then
120			write(default_error,*) 'charge-charge does not have cutoffs set!'
121			return
122			endif
123
124			if (.not.haveChargeMap) then
125			localError = 0
126			call createChargeMap(localError)
127			if ( localError .ne. 0 ) then
128			call handleError("charge-charge", "ChargeMap creation failed!")
129			return
130			end if
131			endif
132
133			#ifdef IS_MPI
134			me1 = atid_Row(atom1)
135			me2 = atid_Col(atom2)
136			#else
137			me1 = atid(atom1)
138			me2 = atid(atom2)
139			#endif
140
141			q1 = ChargeMap(me1)%charge
142			q2 = ChargeMap(me2)%charge
143	gezelter	1138
144	gezelter	1150	vterm = sw * pre * q1 * q2 / rij
145	gezelter	1138
146	gezelter	1150	dudr = -vterm / rij
147			drdx = d(1) / rij
148			drdy = d(2) / rij
149			drdz = d(3) / rij
150
151			fx = dudr * drdx
152			fy = dudr * drdy
153			fz = dudr * drdz
154	gezelter	940
155	gezelter	1150	vpair = vpair + vterm
156
157	gezelter	940	#ifdef IS_MPI
158	gezelter	1150	if (do_pot) then
159			pot_Row(atom1) = pot_Row(atom1) + vterm*0.5
160			pot_Col(atom2) = pot_Col(atom2) + vterm*0.5
161			endif
162
163			f_Row(1,atom1) = f_Row(1,atom1) + fx
164			f_Row(2,atom1) = f_Row(2,atom1) + fy
165			f_Row(3,atom1) = f_Row(3,atom1) + fz
166
167			f_Col(1,atom2) = f_Col(1,atom2) - fx
168			f_Col(2,atom2) = f_Col(2,atom2) - fy
169			f_Col(3,atom2) = f_Col(3,atom2) - fz
170
171	gezelter	940	#else
172
173	gezelter	1150	if (do_pot) pot = pot + vterm
174
175			f(1,atom1) = f(1,atom1) + fx
176			f(2,atom1) = f(2,atom1) + fy
177			f(3,atom1) = f(3,atom1) + fz
178
179			f(1,atom2) = f(1,atom2) - fx
180			f(2,atom2) = f(2,atom2) - fy
181			f(3,atom2) = f(3,atom2) - fz
182
183	gezelter	940	#endif
184	gezelter	1150
185			if (do_stress) then
186
187	gezelter	940	#ifdef IS_MPI
188	gezelter	1150	id1 = tagRow(atom1)
189			id2 = tagColumn(atom2)
190	gezelter	940	#else
191	gezelter	1150	id1 = atom1
192			id2 = atom2
193	gezelter	940	#endif
194	gezelter	1150
195
196			if (molMembershipList(id1) .ne. molMembershipList(id2)) then
197
198			! because the d vector is the rj - ri vector, and
199			! because fx, fy, fz are the force on atom i, we need a
200			! negative sign here:
201
202			tau_Temp(1) = tau_Temp(1) - d(1) * fx
203			tau_Temp(2) = tau_Temp(2) - d(1) * fy
204			tau_Temp(3) = tau_Temp(3) - d(1) * fz
205			tau_Temp(4) = tau_Temp(4) - d(2) * fx
206			tau_Temp(5) = tau_Temp(5) - d(2) * fy
207			tau_Temp(6) = tau_Temp(6) - d(2) * fz
208			tau_Temp(8) = tau_Temp(8) - d(3) * fy
209			tau_Temp(9) = tau_Temp(9) - d(3) * fz
210
211
212			virial_Temp = virial_Temp + &
213			(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
214
215	gezelter	940	endif
216			endif
217	gezelter	1150
218
219	gezelter	940	return
220			end subroutine do_charge_pair
221
222			end module charge_charge