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

Diff Legend

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