ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/calc_charge_charge.F90
Revision: 1129
Committed: Thu Apr 22 03:29:30 2004 UTC (20 years, 2 months ago) by tim
File size: 5123 byte(s)
Log Message:
fixed two bugs in MPI version of InitfromFile and one unmatch MPI command in DumpWriter

File Contents

# 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     real(kind=dp), save :: pre = 0.0_DP
18     logical, save :: haveCutoffs = .false.
19     logical, save :: haveChargeMap = .false.
20    
21 gezelter 945 public::setCutoffsCharge
22 gezelter 940 public::do_charge_pair
23    
24     type :: ChargeList
25     real(kind=DP) :: charge = 0.0_DP
26     end type ChargeList
27    
28     type(ChargeList), dimension(:), allocatable :: ChargeMap
29    
30     contains
31    
32 gezelter 945 subroutine setCutoffsCharge(this_ecr, this_rt)
33     real(kind=dp), intent(in) :: this_ecr, this_rt
34     ecr = this_ecr
35 gezelter 940 rt = this_rt
36    
37     ! pre converts from fundamental charge to kcal/mol
38     pre = 332.06508_DP
39    
40     haveCutoffs = .true.
41    
42     return
43 gezelter 945 end subroutine setCutoffsCharge
44 gezelter 940
45     subroutine createChargeMap(status)
46     integer :: nAtypes
47     integer :: status
48     integer :: i
49     real (kind=DP) :: thisCharge
50     logical :: thisProperty
51    
52     status = 0
53    
54     nAtypes = getSize(atypes)
55    
56     if (nAtypes == 0) then
57     status = -1
58     return
59     end if
60    
61     if (.not. allocated(ChargeMap)) then
62     allocate(ChargeMap(nAtypes))
63     endif
64    
65     do i = 1, nAtypes
66    
67     call getElementProperty(atypes, i, "is_Charge", thisProperty)
68    
69     if (thisProperty) then
70     call getElementProperty(atypes, i, "charge", thisCharge)
71     ChargeMap(i)%charge = thisCharge
72     endif
73    
74     end do
75    
76     haveChargeMap = .true.
77    
78     end subroutine createChargeMap
79    
80    
81     subroutine do_charge_pair(atom1, atom2, d, rij, r2, pot, f, &
82     do_pot, do_stress)
83    
84     logical :: do_pot, do_stress
85    
86     integer atom1, atom2, me1, me2, id1, id2
87     integer :: localError
88     real(kind=dp) :: rij, r2, q1, q2
89     real(kind=dp) :: drdx, drdy, drdz, dudr, fx, fy, fz
90     real(kind=dp) :: taper, dtdr, vterm
91    
92     real( kind = dp ) :: pot
93     real( kind = dp ), dimension(3) :: d
94     real( kind = dp ), dimension(3,nLocal) :: f
95    
96     if (.not. haveCutoffs) then
97     write(default_error,*) 'charge-charge does not have cutoffs set!'
98     return
99     endif
100    
101     if (.not.haveChargeMap) then
102     localError = 0
103     call createChargeMap(localError)
104     if ( localError .ne. 0 ) then
105     call handleError("charge-charge", "ChargeMap creation failed!")
106     return
107     end if
108     endif
109    
110     #ifdef IS_MPI
111     me1 = atid_Row(atom1)
112     me2 = atid_Col(atom2)
113     #else
114     me1 = atid(atom1)
115     me2 = atid(atom2)
116     #endif
117    
118     q1 = ChargeMap(me1)%charge
119     q2 = ChargeMap(me2)%charge
120    
121 gezelter 945 if (rij.le.ecr) then
122 gezelter 940
123     if (rij.lt.rt) then
124     taper = 1.0d0
125     dtdr = 0.0d0
126     else
127 gezelter 945 taper = (ecr + 2.0d0*rij - 3.0d0*rt)*(ecr-rij)**2/ ((ecr-rt)**3)
128     dtdr = 6.0d0*(rij*rij - rij*rt - rij*ecr +ecr*rt)/((ecr-rt)**3)
129 gezelter 940 endif
130    
131     vterm = pre * q1 * q2 / rij
132    
133 tim 1129 !!if(rij .le. 1.5) then
134     !! write(*,*) atom1, atom2, q1, q2, rij, vterm
135     !!endif
136    
137 gezelter 940 dudr = vterm * ( dtdr - taper / rij )
138    
139     drdx = d(1) / rij
140     drdy = d(2) / rij
141     drdz = d(3) / rij
142    
143     fx = dudr * drdx
144     fy = dudr * drdy
145 tim 1113 fz = dudr * drdz
146    
147 gezelter 940 #ifdef IS_MPI
148     if (do_pot) then
149     pot_Row(atom1) = pot_Row(atom1) + vterm*taper*0.5
150     pot_Col(atom2) = pot_Col(atom2) + vterm*taper*0.5
151     endif
152    
153     f_Row(1,atom1) = f_Row(1,atom1) + fx
154     f_Row(2,atom1) = f_Row(2,atom1) + fy
155     f_Row(3,atom1) = f_Row(3,atom1) + fz
156    
157     f_Col(1,atom2) = f_Col(1,atom2) - fx
158     f_Col(2,atom2) = f_Col(2,atom2) - fy
159     f_Col(3,atom2) = f_Col(3,atom2) - fz
160    
161     #else
162    
163     if (do_pot) pot = pot + vterm*taper
164    
165     f(1,atom1) = f(1,atom1) + fx
166     f(2,atom1) = f(2,atom1) + fy
167     f(3,atom1) = f(3,atom1) + fz
168    
169     f(1,atom2) = f(1,atom2) - fx
170     f(2,atom2) = f(2,atom2) - fy
171     f(3,atom2) = f(3,atom2) - fz
172    
173     #endif
174    
175     if (do_stress) then
176    
177     #ifdef IS_MPI
178     id1 = tagRow(atom1)
179     id2 = tagColumn(atom2)
180     #else
181     id1 = atom1
182     id2 = atom2
183     #endif
184    
185     if (molMembershipList(id1) .ne. molMembershipList(id2)) then
186    
187     ! because the d vector is the rj - ri vector, and
188     ! because fx, fy, fz are the force on atom i, we need a
189     ! negative sign here:
190    
191     tau_Temp(1) = tau_Temp(1) - d(1) * fx
192     tau_Temp(2) = tau_Temp(2) - d(1) * fy
193     tau_Temp(3) = tau_Temp(3) - d(1) * fz
194     tau_Temp(4) = tau_Temp(4) - d(2) * fx
195     tau_Temp(5) = tau_Temp(5) - d(2) * fy
196     tau_Temp(6) = tau_Temp(6) - d(2) * fz
197     tau_Temp(7) = tau_Temp(7) - d(3) * fx
198     tau_Temp(8) = tau_Temp(8) - d(3) * fy
199     tau_Temp(9) = tau_Temp(9) - d(3) * fz
200    
201     virial_Temp = virial_Temp + &
202     (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
203    
204     endif
205     endif
206    
207     endif
208     return
209     end subroutine do_charge_pair
210    
211     end module charge_charge