ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-4/src/UseTheForce/doForces.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 1610 by gezelter, Wed Oct 20 04:19:55 2004 UTC vs.
Revision 1634 by gezelter, Fri Oct 22 21:21:02 2004 UTC

# Line 4 | Line 4
4  
5   !! @author Charles F. Vardeman II
6   !! @author Matthew Meineke
7 < !! @version $Id: doForces.F90,v 1.1 2004-10-20 04:19:55 gezelter Exp $, $Date: 2004-10-20 04:19:55 $, $Name: not supported by cvs2svn $, $Revision: 1.1 $
7 > !! @version $Id: doForces.F90,v 1.3 2004-10-22 21:20:53 gezelter Exp $, $Date: 2004-10-22 21:20:53 $, $Name: not supported by cvs2svn $, $Revision: 1.3 $
8  
9   module doForces
10    use force_globals
# Line 30 | Line 30 | module doForces
30    PRIVATE
31  
32   #define __FORTRAN90
33 – #include "UseTheForce/fForceField.h"
33   #include "UseTheForce/fSwitchingFunction.h"
34  
35    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
# Line 38 | Line 37 | module doForces
37  
38    logical, save :: haveRlist = .false.
39    logical, save :: haveNeighborList = .false.
41 –  logical, save :: havePolicies = .false.
40    logical, save :: haveSIMvariables = .false.
41    logical, save :: havePropertyMap = .false.
42    logical, save :: haveSaneForceField = .false.
43 <  logical, save :: FF_uses_LJ
44 <  logical, save :: FF_uses_sticky
43 >  
44 >  logical, save :: FF_uses_DirectionalAtoms
45 >  logical, save :: FF_uses_LennardJones
46 >  logical, save :: FF_uses_Electrostatic
47    logical, save :: FF_uses_charges
48    logical, save :: FF_uses_dipoles
49 <  logical, save :: FF_uses_RF
50 <  logical, save :: FF_uses_GB
49 >  logical, save :: FF_uses_sticky
50 >  logical, save :: FF_uses_GayBerne
51    logical, save :: FF_uses_EAM
52 <  logical, save :: SIM_uses_LJ
53 <  logical, save :: SIM_uses_sticky
54 <  logical, save :: SIM_uses_charges
55 <  logical, save :: SIM_uses_dipoles
56 <  logical, save :: SIM_uses_RF
57 <  logical, save :: SIM_uses_GB
52 >  logical, save :: FF_uses_Shapes
53 >  logical, save :: FF_uses_FLARB
54 >  logical, save :: FF_uses_RF
55 >
56 >  logical, save :: SIM_uses_DirectionalAtoms
57 >  logical, save :: SIM_uses_LennardJones
58 >  logical, save :: SIM_uses_Electrostatics
59 >  logical, save :: SIM_uses_Charges
60 >  logical, save :: SIM_uses_Dipoles
61 >  logical, save :: SIM_uses_Sticky
62 >  logical, save :: SIM_uses_GayBerne
63    logical, save :: SIM_uses_EAM
64 +  logical, save :: SIM_uses_Shapes
65 +  logical, save :: SIM_uses_FLARB
66 +  logical, save :: SIM_uses_RF
67    logical, save :: SIM_requires_postpair_calc
68    logical, save :: SIM_requires_prepair_calc
61 –  logical, save :: SIM_uses_directional_atoms
69    logical, save :: SIM_uses_PBC
70    logical, save :: SIM_uses_molecular_cutoffs
71  
# Line 76 | Line 83 | module doForces
83   #endif
84  
85    type :: Properties
86 <     logical :: is_lj     = .false.
87 <     logical :: is_sticky = .false.
88 <     logical :: is_dp     = .false.
89 <     logical :: is_gb     = .false.
90 <     logical :: is_eam    = .false.
91 <     logical :: is_charge = .false.
92 <     real(kind=DP) :: charge = 0.0_DP
93 <     real(kind=DP) :: dipole_moment = 0.0_DP
86 >     logical :: is_Directional   = .false.
87 >     logical :: is_LennardJones  = .false.
88 >     logical :: is_Electrostatic = .false.
89 >     logical :: is_Charge        = .false.
90 >     logical :: is_Dipole        = .false.
91 >     logical :: is_Sticky        = .false.
92 >     logical :: is_GayBerne      = .false.
93 >     logical :: is_EAM           = .false.
94 >     logical :: is_Shape         = .false.
95 >     logical :: is_FLARB         = .false.
96    end type Properties
97  
98    type(Properties), dimension(:),allocatable :: PropertyMap
# Line 122 | Line 131 | contains
131      endif
132  
133      do i = 1, nAtypes
134 <       call getElementProperty(atypes, i, "is_LJ", thisProperty)
135 <       PropertyMap(i)%is_LJ = thisProperty
134 >       call getElementProperty(atypes, i, "is_Directional", thisProperty)
135 >       PropertyMap(i)%is_Directional = thisProperty
136  
137 +       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
138 +       PropertyMap(i)%is_LennardJones = thisProperty
139 +      
140 +       call getElementProperty(atypes, i, "is_Electrostatic", thisProperty)
141 +       PropertyMap(i)%is_Electrostatic = thisProperty
142 +
143         call getElementProperty(atypes, i, "is_Charge", thisProperty)
144         PropertyMap(i)%is_Charge = thisProperty
145        
146 <       if (thisProperty) then
147 <          call getElementProperty(atypes, i, "charge", thisDPproperty)
133 <          PropertyMap(i)%charge = thisDPproperty
134 <       endif
135 <
136 <       call getElementProperty(atypes, i, "is_DP", thisProperty)
137 <       PropertyMap(i)%is_DP = thisProperty
138 <
139 <       if (thisProperty) then
140 <          call getElementProperty(atypes, i, "dipole_moment", thisDPproperty)
141 <          PropertyMap(i)%dipole_moment = thisDPproperty
142 <       endif
146 >       call getElementProperty(atypes, i, "is_Dipole", thisProperty)
147 >       PropertyMap(i)%is_Dipole = thisProperty
148  
149         call getElementProperty(atypes, i, "is_Sticky", thisProperty)
150         PropertyMap(i)%is_Sticky = thisProperty
151 <       call getElementProperty(atypes, i, "is_GB", thisProperty)
152 <       PropertyMap(i)%is_GB = thisProperty
151 >
152 >       call getElementProperty(atypes, i, "is_GayBerne", thisProperty)
153 >       PropertyMap(i)%is_GayBerne = thisProperty
154 >
155         call getElementProperty(atypes, i, "is_EAM", thisProperty)
156         PropertyMap(i)%is_EAM = thisProperty
157 +
158 +       call getElementProperty(atypes, i, "is_Shape", thisProperty)
159 +       PropertyMap(i)%is_Shape = thisProperty
160 +
161 +       call getElementProperty(atypes, i, "is_FLARB", thisProperty)
162 +       PropertyMap(i)%is_FLARB = thisProperty
163      end do
164  
165      havePropertyMap = .true.
# Line 154 | Line 167 | contains
167    end subroutine createPropertyMap
168  
169    subroutine setSimVariables()
170 <    SIM_uses_LJ = SimUsesLJ()
171 <    SIM_uses_sticky = SimUsesSticky()
172 <    SIM_uses_charges = SimUsesCharges()
173 <    SIM_uses_dipoles = SimUsesDipoles()
174 <    SIM_uses_RF = SimUsesRF()
175 <    SIM_uses_GB = SimUsesGB()
170 >    SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
171 >    SIM_uses_LennardJones = SimUsesLennardJones()
172 >    SIM_uses_Electrostatics = SimUsesElectrostatics()
173 >    SIM_uses_Charges = SimUsesCharges()
174 >    SIM_uses_Dipoles = SimUsesDipoles()
175 >    SIM_uses_Sticky = SimUsesSticky()
176 >    SIM_uses_GayBerne = SimUsesGayBerne()
177      SIM_uses_EAM = SimUsesEAM()
178 +    SIM_uses_Shapes = SimUsesShapes()
179 +    SIM_uses_FLARB = SimUsesFLARB()
180 +    SIM_uses_RF = SimUsesRF()
181      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
182      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
166 –    SIM_uses_directional_atoms = SimUsesDirectionalAtoms()
183      SIM_uses_PBC = SimUsesPBC()
168 –    !SIM_uses_molecular_cutoffs = SimUsesMolecularCutoffs()
184  
185      haveSIMvariables = .true.
186  
# Line 202 | Line 217 | contains
217         return
218      endif
219  
205 –    if (SIM_uses_LJ .and. FF_uses_LJ) then
206 –       if (.not. havePolicies) then
207 –          write(default_error, *) 'LJ mixing Policies have not been set in doForces!'
208 –          error = -1
209 –          return
210 –       endif
211 –    endif
212 –
220      if (.not. haveNeighborList) then
221         write(default_error, *) 'neighbor list has not been initialized in doForces!'
222         error = -1
# Line 233 | Line 240 | contains
240    end subroutine doReadyCheck
241      
242  
243 <  subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
243 >  subroutine init_FF(use_RF_c, thisStat)
244  
238 –    integer, intent(in) :: LJMIXPOLICY
245      logical, intent(in) :: use_RF_c
246  
247      integer, intent(out) :: thisStat  
# Line 255 | Line 261 | contains
261      !! this will scan through the known atypes and figure out what
262      !! interactions are used by the force field.    
263    
264 <    FF_uses_LJ = .false.
265 <    FF_uses_sticky = .false.
266 <    FF_uses_charges = .false.
267 <    FF_uses_dipoles = .false.
268 <    FF_uses_GB = .false.
264 >    FF_uses_DirectionalAtoms = .false.
265 >    FF_uses_LennardJones = .false.
266 >    FF_uses_Electrostatic = .false.
267 >    FF_uses_Charges = .false.    
268 >    FF_uses_Dipoles = .false.
269 >    FF_uses_Sticky = .false.
270 >    FF_uses_GayBerne = .false.
271      FF_uses_EAM = .false.
272 +    FF_uses_Shapes = .false.
273 +    FF_uses_FLARB = .false.
274      
275 <    call getMatchingElementList(atypes, "is_LJ", .true., nMatches, MatchList)
276 <    if (nMatches .gt. 0) FF_uses_LJ = .true.
275 >    call getMatchingElementList(atypes, "is_Directional", .true., &
276 >         nMatches, MatchList)
277 >    if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
278  
279 <    call getMatchingElementList(atypes, "is_Charge", .true., nMatches, MatchList)
280 <    if (nMatches .gt. 0) FF_uses_charges = .true.  
279 >    call getMatchingElementList(atypes, "is_LennardJones", .true., &
280 >         nMatches, MatchList)
281 >    if (nMatches .gt. 0) FF_uses_LennardJones = .true.
282 >    
283 >    call getMatchingElementList(atypes, "is_Electrostatic", .true., &
284 >         nMatches, MatchList)
285 >    if (nMatches .gt. 0) then
286 >       FF_uses_Electrostatic = .true.
287 >    endif
288  
289 <    call getMatchingElementList(atypes, "is_DP", .true., nMatches, MatchList)
290 <    if (nMatches .gt. 0) FF_uses_dipoles = .true.
289 >    call getMatchingElementList(atypes, "is_Charge", .true., &
290 >         nMatches, MatchList)
291 >    if (nMatches .gt. 0) then
292 >       FF_uses_charges = .true.  
293 >       FF_uses_electrostatic = .true.
294 >    endif
295      
296 +    call getMatchingElementList(atypes, "is_Dipole", .true., &
297 +         nMatches, MatchList)
298 +    if (nMatches .gt. 0) then
299 +       FF_uses_dipoles = .true.
300 +       FF_uses_electrostatic = .true.
301 +       FF_uses_DirectionalAtoms = .true.
302 +    endif
303 +    
304      call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
305           MatchList)
306 <    if (nMatches .gt. 0) FF_uses_Sticky = .true.
306 >    if (nMatches .gt. 0) then
307 >       FF_uses_Sticky = .true.
308 >       FF_uses_DirectionalAtoms = .true.
309 >    endif
310      
311 <    call getMatchingElementList(atypes, "is_GB", .true., nMatches, MatchList)
312 <    if (nMatches .gt. 0) FF_uses_GB = .true.
311 >    call getMatchingElementList(atypes, "is_GayBerne", .true., &
312 >         nMatches, MatchList)
313 >    if (nMatches .gt. 0) then
314 >       FF_uses_GayBerne = .true.
315 >       FF_uses_DirectionalAtoms = .true.
316 >    endif
317      
318      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
319      if (nMatches .gt. 0) FF_uses_EAM = .true.
320      
321 +    call getMatchingElementList(atypes, "is_Shape", .true., &
322 +         nMatches, MatchList)
323 +    if (nMatches .gt. 0) then
324 +       FF_uses_Shapes = .true.
325 +       FF_uses_DirectionalAtoms = .true.
326 +    endif
327 +
328 +    call getMatchingElementList(atypes, "is_FLARB", .true., &
329 +         nMatches, MatchList)
330 +    if (nMatches .gt. 0) FF_uses_FLARB = .true.
331 +
332      !! Assume sanity (for the sake of argument)
333      haveSaneForceField = .true.
334 <
334 >    
335      !! check to make sure the FF_uses_RF setting makes sense
336      
337      if (FF_uses_dipoles) then
# Line 300 | Line 348 | contains
348         endif
349      endif
350  
303 –    if (FF_uses_LJ) then
304 –      
305 –       select case (LJMIXPOLICY)
306 –       case (LB_MIXING_RULE)
307 –          call init_lj_FF(LB_MIXING_RULE, my_status)            
308 –       case (EXPLICIT_MIXING_RULE)
309 –          call init_lj_FF(EXPLICIT_MIXING_RULE, my_status)
310 –       case default
311 –          write(default_error,*) 'unknown LJ Mixing Policy!'
312 –          thisStat = -1
313 –          haveSaneForceField = .false.
314 –          return            
315 –       end select
316 –       if (my_status /= 0) then
317 –          thisStat = -1
318 –          haveSaneForceField = .false.
319 –          return
320 –       end if
321 –       havePolicies = .true.
322 –    endif
323 –
351      if (FF_uses_sticky) then
352         call check_sticky_FF(my_status)
353         if (my_status /= 0) then
# Line 330 | Line 357 | contains
357         end if
358      endif
359  
333 –
360      if (FF_uses_EAM) then
361           call init_EAM_FF(my_status)
362         if (my_status /= 0) then
# Line 341 | Line 367 | contains
367         end if
368      endif
369  
370 <    if (FF_uses_GB) then
370 >    if (FF_uses_GayBerne) then
371         call check_gb_pair_FF(my_status)
372         if (my_status .ne. 0) then
373            thisStat = -1
# Line 350 | Line 376 | contains
376         endif
377      endif
378  
379 <    if (FF_uses_GB .and. FF_uses_LJ) then
379 >    if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
380      endif
381 +    
382      if (.not. haveNeighborList) then
383         !! Create neighbor lists
384         call expandNeighborList(nLocal, my_status)
# Line 361 | Line 388 | contains
388            return
389         endif
390         haveNeighborList = .true.
391 <    endif
365 <
366 <    
391 >    endif    
392      
393    end subroutine init_FF
394    
# Line 453 | Line 478 | contains
478      call gather(q_group, q_group_Row, plan_group_row_3d)
479      call gather(q_group, q_group_Col, plan_group_col_3d)
480          
481 <    if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
481 >    if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then
482         call gather(u_l, u_l_Row, plan_atom_row_3d)
483         call gather(u_l, u_l_Col, plan_atom_col_3d)
484        
# Line 707 | Line 732 | contains
732         f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
733      end do
734      
735 <    if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
735 >    if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then
736         t_temp = 0.0_dp
737         call scatter(t_Row,t_temp,plan_atom_row_3d)
738         do i = 1,nlocal
# Line 762 | Line 787 | contains
787               me_i = atid(i)
788   #endif
789              
790 <             if (PropertyMap(me_i)%is_DP) then
790 >             if (PropertyMap(me_i)%is_Dipole) then
791                  
792 <                mu_i = PropertyMap(me_i)%dipole_moment
792 >                mu_i = getDipoleMoment(me_i)
793                  
794                  !! The reaction field needs to include a self contribution
795                  !! to the field:
# Line 840 | Line 865 | contains
865      me_j = atid(j)
866   #endif
867      
868 <    if (FF_uses_LJ .and. SIM_uses_LJ) then
868 >    if (FF_uses_LennardJones .and. SIM_uses_LennardJones) then
869        
870 <       if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) then
870 >       if ( PropertyMap(me_i)%is_LennardJones .and. &
871 >            PropertyMap(me_j)%is_LennardJones ) then
872            call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
873         endif
874        
# Line 851 | Line 877 | contains
877      if (FF_uses_charges .and. SIM_uses_charges) then
878        
879         if (PropertyMap(me_i)%is_Charge .and. PropertyMap(me_j)%is_Charge) then
880 <          call do_charge_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
880 >          call do_charge_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
881 >               pot, f, do_pot)
882         endif
883        
884      endif
885      
886      if (FF_uses_dipoles .and. SIM_uses_dipoles) then
887        
888 <       if ( PropertyMap(me_i)%is_DP .and. PropertyMap(me_j)%is_DP) then
889 <          call do_dipole_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, u_l, f, t, &
890 <               do_pot)
888 >       if ( PropertyMap(me_i)%is_Dipole .and. PropertyMap(me_j)%is_Dipole) then
889 >          call do_dipole_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
890 >               pot, u_l, f, t, do_pot)
891            if (FF_uses_RF .and. SIM_uses_RF) then
892               call accumulate_rf(i, j, r, u_l, sw)
893               call rf_correct_forces(i, j, d, r, u_l, sw, f, fpair)
894 <          endif          
894 >          endif
895         endif
896  
897      endif
# Line 872 | Line 899 | contains
899      if (FF_uses_Sticky .and. SIM_uses_sticky) then
900  
901         if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
902 <          call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, A, f, t, &
903 <               do_pot)
902 >          call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
903 >               pot, A, f, t, do_pot)
904         endif
905 <
905 >      
906      endif
907  
908  
909 <    if (FF_uses_GB .and. SIM_uses_GB) then
910 <      
911 <       if ( PropertyMap(me_i)%is_GB .and. PropertyMap(me_j)%is_GB) then
912 <          call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, u_l, f, t, &
913 <               do_pot)
909 >    if (FF_uses_GayBerne .and. SIM_uses_GayBerne) then
910 >      
911 >       if ( PropertyMap(me_i)%is_GayBerne .and. &
912 >            PropertyMap(me_j)%is_GayBerne) then
913 >          call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
914 >               pot, u_l, f, t, do_pot)
915         endif
888 –
889 –    endif
916        
917 +    endif
918 +    
919      if (FF_uses_EAM .and. SIM_uses_EAM) then
920        
921         if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
922            call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
923                 do_pot)
924 +       endif
925 +      
926 +    endif
927 +
928 +    if (FF_uses_Shapes .and. SIM_uses_Shapes) then
929 +      
930 +       if ( PropertyMap(me_i)%is_Shape .and. &
931 +            PropertyMap(me_j)%is_Shape ) then
932 +          call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
933 +               pot, u_l, f, t, do_pot)
934         endif
935        
936      endif
# Line 1135 | Line 1173 | contains
1173  
1174   function FF_UsesDirectionalAtoms() result(doesit)
1175     logical :: doesit
1176 <   doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
1177 <        FF_uses_GB .or. FF_uses_RF
1176 >   doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1177 >        FF_uses_Sticky .or. FF_uses_GayBerne .or. FF_uses_Shapes
1178   end function FF_UsesDirectionalAtoms
1179  
1180   function FF_RequiresPrepairCalc() result(doesit)
# Line 1185 | Line 1223 | end module doForces
1223  
1224   !! Interfaces for C programs to module....
1225  
1226 < subroutine initFortranFF(LJMIXPOLICY, use_RF_c, thisStat)
1226 > subroutine initFortranFF(use_RF_c, thisStat)
1227      use doForces, ONLY: init_FF
1190 –    integer, intent(in) :: LJMIXPOLICY
1228      logical, intent(in) :: use_RF_c
1229  
1230      integer, intent(out) :: thisStat  
1231 <    call init_FF(LJMIXPOLICY, use_RF_c, thisStat)
1231 >    call init_FF(use_RF_c, thisStat)
1232  
1233   end subroutine initFortranFF
1234  

Diff Legend

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