ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/calc_reaction_field.F90
(Generate patch)

Comparing trunk/OOPSE/libmdtools/calc_reaction_field.F90 (file contents):
Revision 898 by chuckv, Mon Jan 5 22:49:14 2004 UTC vs.
Revision 901 by chuckv, Tue Jan 6 19:49:18 2004 UTC

# Line 4 | Line 4 | module reaction_field
4    use atype_module
5    use vector_class
6    use simulation
7 +  use status
8   #ifdef IS_MPI
9    use mpiSimulation
10   #endif
# Line 16 | Line 17 | module reaction_field
17    real(kind=dp), save :: dielect = 1.0_dp
18    real(kind=dp), save :: rrfsq = 1.0_dp
19    real(kind=dp), save :: pre
20 <  logical, save :: rf_initialized = .false., haveCuts = .false.
21 <  logical, save :: haveDie = .false.
20 >  logical, save :: haveCutoffs = .false.
21 >  logical, save :: haveMomentMap = .false.
22 >  logical, save :: haveDielectric = .false.
23  
24 +  type :: MomentList
25 +     real(kind=DP) :: dipole_moment = 0.0_DP
26 +  end type MomentList
27 +
28 +  type(MomentList), dimension(:),allocatable :: MomentMap
29 +
30    PUBLIC::initialize_rf
31    PUBLIC::setCutoffsRF
32    PUBLIC::accumulate_rf
# Line 35 | Line 43 | contains
43  
44      pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf)
45      
46 <    haveDie = .true.
39 <    if (haveCuts) rf_initialized = .true.
46 >    haveDielectric = .true.
47  
48      return
49    end subroutine initialize_rf
# Line 51 | Line 58 | contains
58      rrfsq = rrf * rrf
59      pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf)
60      
61 <    haveCuts = .true.
55 <    if (haveDie) rf_initialized = .true.
61 >    haveCutoffs = .true.
62  
63    end subroutine setCutoffsRF
64  
65 <  
65 >  subroutine createMomentMap(status)
66 >    integer :: nAtypes
67 >    integer :: status
68 >    integer :: i
69 >    real (kind=DP) :: thisDP
70 >    logical :: thisProperty
71 >
72 >    status = 0
73 >
74 >    nAtypes = getSize(atypes)
75 >    
76 >    if (nAtypes == 0) then
77 >       status = -1
78 >       return
79 >    end if
80 >    
81 >    if (.not. allocated(MomentMap)) then
82 >       allocate(MomentMap(nAtypes))
83 >    endif
84 >
85 >    do i = 1, nAtypes
86 >
87 >       call getElementProperty(atypes, i, "is_DP", thisProperty)
88 >
89 >       if (thisProperty) then
90 >          call getElementProperty(atypes, i, "dipole_moment", thisDP)
91 >          MomentMap(i)%dipole_moment = thisDP
92 >       endif
93 >      
94 >    end do
95 >    
96 >    haveMomentMap = .true.
97 >    
98 >  end subroutine createMomentMap  
99 >
100    subroutine accumulate_rf(atom1, atom2, rij, u_l)
101  
102      integer, intent(in) :: atom1, atom2
# Line 68 | Line 108 | contains
108      real (kind = dp), dimension(3) :: ul1
109      real (kind = dp), dimension(3) :: ul2  
110  
111 <    if (.not.rf_initialized) then
111 >    integer :: localError
112 >
113 >    if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
114         write(default_error,*) 'Reaction field not initialized!'
115         return
116      endif
117 +
118 +    if (.not.haveMomentMap) then
119 +       localError = 0
120 +       call createMomentMap(localError)
121 +       if ( localError .ne. 0 ) then
122 +          call handleError("reaction-field", "MomentMap creation failed!")
123 +          return
124 +       end if
125 +    endif
126      
127      if (rij.le.rrf) then
128            
# Line 104 | Line 155 | contains
155         ul2(3) = u_l(3,atom2)
156   #endif
157        
158 <       call getElementProperty(atypes, me1, "dipole_moment", mu1)
159 <       call getElementProperty(atypes, me2, "dipole_moment", mu2)
158 >       mu1 = MomentMap(me1)%dipole_moment
159 >       mu2 = MomentMap(me2)%dipole_moment
160        
110 –      
161   #ifdef IS_MPI
162         rf_Row(1,atom1) = rf_Row(1,atom1) + ul2(1)*mu2*taper
163         rf_Row(2,atom1) = rf_Row(2,atom1) + ul2(2)*mu2*taper
# Line 153 | Line 203 | contains
203      real (kind = dp), dimension(3,nLocal) :: u_l    
204      real (kind = dp), dimension(3,nLocal) :: t
205  
206 <    if (.not.rf_initialized) then
206 >    integer :: localError
207 >
208 >    if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
209         write(default_error,*) 'Reaction field not initialized!'
210         return
211      endif
212  
213 +    if (.not.haveMomentMap) then
214 +       localError = 0
215 +       call createMomentMap(localError)
216 +       if ( localError .ne. 0 ) then
217 +          call handleError("reaction-field", "MomentMap creation failed!")
218 +          return
219 +       end if
220 +    endif
221 +
222      ! compute torques on dipoles:
223      ! pre converts from mu in units of debye to kcal/mol
224      
# Line 193 | Line 254 | contains
254      integer :: me1, me2, id1, id2
255      real (kind = dp) :: mu1, mu2
256      
257 <    if (.not.rf_initialized) then
257 >    integer :: localError
258 >
259 >    if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
260         write(default_error,*) 'Reaction field not initialized!'
261         return
262      endif
263  
264 +    if (.not.haveMomentMap) then
265 +       localError = 0
266 +       call createMomentMap(localError)
267 +       if ( localError .ne. 0 ) then
268 +          call handleError("reaction-field", "MomentMap creation failed!")
269 +          return
270 +       end if
271 +    endif
272 +
273      if (rij.le.rrf) then
274        
275         if (rij.lt.rt) then
# Line 229 | Line 301 | contains
301         ul2(3) = u_l(3,atom2)
302   #endif
303        
304 <       call getElementProperty(atypes, me1, "dipole_moment", mu1)
305 <       call getElementProperty(atypes, me2, "dipole_moment", mu2)
304 >       mu1 = MomentMap(me1)%dipole_moment
305 >       mu2 = MomentMap(me2)%dipole_moment
306        
307         u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)
308        

Diff Legend

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