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

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

# Line 4 | Line 4
4  
5   !! @author Charles F. Vardeman II
6   !! @author Matthew Meineke
7 < !! @version $Id: do_Forces.F90,v 1.44 2004-01-05 22:49:14 chuckv Exp $, $Date: 2004-01-05 22:49:14 $, $Name: not supported by cvs2svn $, $Revision: 1.44 $
7 > !! @version $Id: do_Forces.F90,v 1.45 2004-01-06 18:54:57 chuckv Exp $, $Date: 2004-01-06 18:54:57 $, $Name: not supported by cvs2svn $, $Revision: 1.45 $
8  
9   module do_Forces
10    use force_globals
# Line 30 | Line 30 | module do_Forces
30   #define __FORTRAN90
31   #include "fForceField.h"
32  
33 <  logical, save :: do_forces_initialized = .false., haveRlist = .false.
33 >  logical, save :: haveRlist = .false.
34 >  logical, save :: haveNeighborList = .false.
35    logical, save :: havePolicies = .false.
36 +  logical, save :: haveSIMvariables = .false.
37 +  logical, save :: havePropertyMap = .false.
38 +  logical, save :: haveSaneForceField = .false.
39    logical, save :: FF_uses_LJ
40    logical, save :: FF_uses_sticky
41    logical, save :: FF_uses_dipoles
42    logical, save :: FF_uses_RF
43    logical, save :: FF_uses_GB
44    logical, save :: FF_uses_EAM
45 +  logical, save :: SIM_uses_LJ
46 +  logical, save :: SIM_uses_sticky
47 +  logical, save :: SIM_uses_dipoles
48 +  logical, save :: SIM_uses_RF
49 +  logical, save :: SIM_uses_GB
50 +  logical, save :: SIM_uses_EAM
51 +  logical, save :: SIM_requires_postpair_calc
52 +  logical, save :: SIM_requires_prepair_calc
53 +  logical, save :: SIM_uses_directional_atoms
54 +  logical, save :: SIM_uses_PBC
55  
56    real(kind=dp), save :: rlist, rlistsq
57  
# Line 45 | Line 59 | module do_Forces
59    public :: do_force_loop
60    public :: setRlistDF
61  
62 +
63   #ifdef PROFILE
64    public :: getforcetime
65    real, save :: forceTime = 0
# Line 52 | Line 67 | module do_Forces
67    integer :: nLoops
68   #endif
69  
70 <  logical, allocatable :: propertyMapI(:,:)
71 <  logical, allocatable :: propertyMapJ(:,:)
70 >  type :: Properties
71 >     logical :: is_lj     = .false.
72 >     logical :: is_sticky = .false.
73 >     logical :: is_dp     = .false.
74 >     logical :: is_gb     = .false.
75 >     logical :: is_eam    = .false.
76 >     real(kind=DP) :: dipole_moment = 0.0_DP
77 >  end type Properties
78  
79 +  type(Properties), dimension(:),allocatable :: PropertyMap
80 +
81   contains
82  
83    subroutine setRlistDF( this_rlist )
# Line 65 | Line 88 | contains
88      rlistsq = rlist * rlist
89      
90      haveRlist = .true.
68    if( havePolicies ) do_forces_initialized = .true.
91  
92    end subroutine setRlistDF    
93 +
94 +  subroutine createPropertyMap(status)
95 +    integer :: nAtypes
96 +    integer :: status
97 +    integer :: i
98 +    logical :: thisProperty
99 +    real (kind=DP) :: thisDPproperty
100 +
101 +    status = 0
102 +
103 +    nAtypes = getSize(atypes)
104 +
105 +    if (nAtypes == 0) then
106 +       status = -1
107 +       return
108 +    end if
109 +        
110 +    if (.not. allocated(PropertyMap)) then
111 +       allocate(PropertyMap(nAtypes))
112 +    endif
113 +
114 +    do i = 1, nAtypes
115 +       call getElementProperty(atypes, i, "is_LJ", thisProperty)
116 +       PropertyMap(i)%is_LJ = thisProperty
117 +       call getElementProperty(atypes, i, "is_DP", thisProperty)
118 +       PropertyMap(i)%is_DP = thisProperty
119 +
120 +       if (thisProperty) then
121 +          call getElementProperty(atypes, i, "dipole_moment", thisDPproperty)
122 +          PropertyMap(i)%dipole_moment = thisDPproperty
123 +       endif
124 +
125 +       call getElementProperty(atypes, i, "is_Sticky", thisProperty)
126 +       PropertyMap(i)%is_Sticky = thisProperty
127 +       call getElementProperty(atypes, i, "is_GB", thisProperty)
128 +       PropertyMap(i)%is_GB = thisProperty
129 +       call getElementProperty(atypes, i, "is_EAM", thisProperty)
130 +       PropertyMap(i)%is_EAM = thisProperty
131 +    end do
132 +
133 +    havePropertyMap = .true.
134 +
135 +  end subroutine createPropertyMap
136 +
137 +  subroutine setSimVariables()
138 +    SIM_uses_LJ = SimUsesLJ()
139 +    SIM_uses_sticky = SimUsesSticky()
140 +    SIM_uses_dipoles = SimUsesDipoles()
141 +    SIM_uses_RF = SimUsesRF()
142 +    SIM_uses_GB = SimUsesGB()
143 +    SIM_uses_EAM = SimUsesEAM()
144 +    SIM_requires_postpair_calc = SimRequiresPostpairCalc()
145 +    SIM_requires_prepair_calc = SimRequiresPrepairCalc()
146 +    SIM_uses_directional_atoms = SimUsesDirectionalAtoms()
147 +    SIM_uses_PBC = SimUsesPBC()
148 +
149 +    haveSIMvariables = .true.
150  
151 +    return
152 +  end subroutine setSimVariables
153 +
154 +  subroutine doReadyCheck(error)
155 +    integer, intent(out) :: error
156 +
157 +    integer :: myStatus
158 +
159 +    error = 0
160 +    
161 +    if (.not. havePropertyMap) then
162 +
163 +       myStatus = 0
164 +
165 +       call createPropertyMap(myStatus)
166 +
167 +       if (myStatus .ne. 0) then
168 +          write(default_error, *) 'createPropertyMap failed in do_Forces!'
169 +          error = -1
170 +          return
171 +       endif
172 +    endif
173 +
174 +    if (.not. haveSIMvariables) then
175 +       call setSimVariables()
176 +    endif
177 +
178 +    if (.not. haveRlist) then
179 +       write(default_error, *) 'rList has not been set in do_Forces!'
180 +       error = -1
181 +       return
182 +    endif
183 +
184 +    if (SIM_uses_LJ .and. FF_uses_LJ) then
185 +       if (.not. havePolicies) then
186 +          write(default_error, *) 'LJ mixing Policies have not been set in do_Forces!'
187 +          error = -1
188 +          return
189 +       endif
190 +    endif
191 +
192 +    if (.not. haveNeighborList) then
193 +       write(default_error, *) 'neighbor list has not been initialized in do_Forces!'
194 +       error = -1
195 +       return
196 +    end if
197 +
198 +    if (.not. haveSaneForceField) then
199 +       write(default_error, *) 'Force Field is not sane in do_Forces!'
200 +       error = -1
201 +       return
202 +    end if
203 +
204 + #ifdef IS_MPI
205 +    if (.not. isMPISimSet()) then
206 +       write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
207 +       error = -1
208 +       return
209 +    endif
210 + #endif
211 +    return
212 +  end subroutine doReadyCheck
213 +    
214 +
215    subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
216  
217      integer, intent(in) :: LJMIXPOLICY
# Line 113 | Line 256 | contains
256      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
257      if (nMatches .gt. 0) FF_uses_EAM = .true.
258      
259 +    !! Assume sanity (for the sake of argument)
260 +    haveSaneForceField = .true.
261 +
262      !! check to make sure the FF_uses_RF setting makes sense
263      
264      if (FF_uses_dipoles) then
# Line 124 | Line 270 | contains
270         if (FF_uses_RF) then          
271            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
272            thisStat = -1
273 +          haveSaneForceField = .false.
274            return
275         endif
276      endif
# Line 138 | Line 285 | contains
285         case default
286            write(default_error,*) 'unknown LJ Mixing Policy!'
287            thisStat = -1
288 +          haveSaneForceField = .false.
289            return            
290         end select
291         if (my_status /= 0) then
292            thisStat = -1
293 +          haveSaneForceField = .false.
294            return
295         end if
296 +       havePolicies = .true.
297      endif
298  
299      if (FF_uses_sticky) then
300         call check_sticky_FF(my_status)
301         if (my_status /= 0) then
302            thisStat = -1
303 +          haveSaneForceField = .false.
304            return
305         end if
306      endif
# Line 158 | Line 309 | contains
309      if (FF_uses_EAM) then
310           call init_EAM_FF(my_status)
311         if (my_status /= 0) then
312 <          write(*,*) "init_EAM_FF returned a bad status"
312 >          write(default_error, *) "init_EAM_FF returned a bad status"
313            thisStat = -1
314 +          haveSaneForceField = .false.
315            return
316         end if
317      endif
318  
167
168    
319      if (FF_uses_GB) then
320         call check_gb_pair_FF(my_status)
321         if (my_status .ne. 0) then
322            thisStat = -1
323 +          haveSaneForceField = .false.
324            return
325         endif
326      endif
327  
328      if (FF_uses_GB .and. FF_uses_LJ) then
329      endif
330 <    if (.not. do_forces_initialized) then
330 >    if (.not. haveNeighborList) then
331         !! Create neighbor lists
332         call expandNeighborList(nLocal, my_status)
333         if (my_Status /= 0) then
# Line 184 | Line 335 | contains
335            thisStat = -1
336            return
337         endif
338 +       haveNeighborList = .true.
339      endif
340      
189
190    havePolicies = .true.
191    if( haveRlist ) do_forces_initialized = .true.
192
341    end subroutine init_FF
342    
343  
# Line 246 | Line 394 | contains
394      natoms = nlocal
395   #endif
396  
397 <    call check_initialization(localError)
397 >    call doReadyCheck(localError)
398      if ( localError .ne. 0 ) then
399 <       call handleError("do_force_loop","Not Initialized")
399 >       call handleError("do_force_loop", "Not Initialized")
400         error = -1
401         return
402      end if
# Line 256 | Line 404 | contains
404  
405      do_pot = do_pot_c
406      do_stress = do_stress_c
259
260
261 #ifdef IS_MPI
262    if (.not.allocated(propertyMapI)) then
263       allocate(propertyMapI(5,nrow))
264    endif
265
266    do i = 1, nrow
267       me_i = atid_row(i)
268 #else
269    if (.not.allocated(propertyMapI)) then
270       allocate(propertyMapI(5,nlocal))
271    endif
272
273    do i = 1, natoms
274       me_i = atid(i)
275 #endif
276      
277       propertyMapI(1:5,i) = .false.
278
279       call getElementProperty(atypes, me_i, "propertyPack", propPack_i)
280    
281       ! unpack the properties
282      
283       if (iand(propPack_i, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) &
284            propertyMapI(1, i) = .true.
285       if (iand(propPack_i, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) &
286            propertyMapI(2, i) = .true.
287       if (iand(propPack_i, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) &
288            propertyMapI(3, i) = .true.
289       if (iand(propPack_i, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) &
290            propertyMapI(4, i) = .true.
291       if (iand(propPack_i, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) &
292            propertyMapI(5, i) = .true.
293
294    end do
295
296 #ifdef IS_MPI
297    if (.not.allocated(propertyMapJ)) then
298       allocate(propertyMapJ(5,ncol))
299    endif
407  
301    do j = 1, ncol
302       me_j = atid_col(j)
303 #else
304    if (.not.allocated(propertyMapJ)) then
305       allocate(propertyMapJ(5,nlocal))
306    endif
307
308    do j = 1, natoms
309       me_j = atid(j)
310 #endif
311      
312       propertyMapJ(1:5,j) = .false.
313
314       call getElementProperty(atypes, me_j, "propertyPack", propPack_j)
315    
316       ! unpack the properties
317      
318       if (iand(propPack_j, LJ_PROPERTY_MASK) .eq. LJ_PROPERTY_MASK) &
319            propertyMapJ(1, j) = .true.
320       if (iand(propPack_j, DP_PROPERTY_MASK) .eq. DP_PROPERTY_MASK) &
321            propertyMapJ(2, j) = .true.
322       if (iand(propPack_j, STICKY_PROPERTY_MASK) .eq. STICKY_PROPERTY_MASK) &
323            propertyMapJ(3, j) = .true.
324       if (iand(propPack_j, GB_PROPERTY_MASK) .eq. GB_PROPERTY_MASK) &
325            propertyMapJ(4, j) = .true.
326       if (iand(propPack_j, EAM_PROPERTY_MASK) .eq. EAM_PROPERTY_MASK) &
327            propertyMapJ(5, j) = .true.
328
329    end do
330
408      ! Gather all information needed by all force loops:
409      
410   #ifdef IS_MPI    
# Line 335 | Line 412 | contains
412      call gather(q,q_Row,plan_row3d)
413      call gather(q,q_Col,plan_col3d)
414          
415 <    if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
415 >    if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
416         call gather(u_l,u_l_Row,plan_row3d)
417         call gather(u_l,u_l_Col,plan_col3d)
418        
# Line 351 | Line 428 | contains
428      nloops = nloops + 1
429   #endif
430    
431 <    if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
431 >    if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
432         !! See if we need to update neighbor lists
433         call checkNeighborList(nlocal, q, listSkin, update_nlist)  
434         !! if_mpi_gather_stuff_for_prepair
# Line 672 | Line 749 | contains
749         f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
750      end do
751      
752 <    if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
752 >    if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
753         t_temp = 0.0_dp
754         call scatter(t_Row,t_temp,plan_row3d)
755         do i = 1,nlocal
# Line 706 | Line 783 | contains
783      endif    
784   #endif
785  
786 <    if (FF_RequiresPostpairCalc() .and. SimRequiresPostpairCalc()) then
786 >    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
787        
788 <       if (FF_uses_RF .and. SimUsesRF()) then
788 >       if (FF_uses_RF .and. SIM_uses_RF) then
789            
790   #ifdef IS_MPI
791            call scatter(rf_Row,rf,plan_row3d)
# Line 726 | Line 803 | contains
803   #else
804               me_i = atid(i)
805   #endif
806 <             call getElementProperty(atypes, me_i, "is_DP", is_DP_i)      
807 <             if ( is_DP_i ) then
808 <                call getElementProperty(atypes, me_i, "dipole_moment", mu_i)
806 >
807 >             if (PropertyMap(me_i)%is_DP) then
808 >
809 >                mu_i = PropertyMap(me_i)%dipole_moment
810 >
811                  !! The reaction field needs to include a self contribution
812                  !! to the field:
813 <                call accumulate_self_rf(i, mu_i, u_l)            
813 >                call accumulate_self_rf(i, mu_i, u_l)
814                  !! Get the reaction field contribution to the
815                  !! potential and torques:
816                  call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot)
# Line 813 | Line 892 | contains
892  
893   #endif
894      
895 <    if (FF_uses_LJ .and. SimUsesLJ()) then
895 >    if (FF_uses_LJ .and. SIM_uses_LJ) then
896  
897 <       if ( propertyMapI(1, me_i) .and. propertyMapJ(1, me_j) ) &
897 >       if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) &
898              call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
899  
900      endif
901  
902 <    if (FF_uses_dipoles .and. SimUsesDipoles()) then
902 >    if (FF_uses_dipoles .and. SIM_uses_dipoles) then
903        
904 <       if ( propertyMapI(2, me_i) .and. propertyMapJ(2, me_j)) then
904 >       if ( PropertyMap(me_i)%is_DP .and. PropertyMap(me_j)%is_DP) then
905            call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
906                 do_pot, do_stress)
907 <          if (FF_uses_RF .and. SimUsesRF()) then
907 >          if (FF_uses_RF .and. SIM_uses_RF) then
908               call accumulate_rf(i, j, r, u_l)
909               call rf_correct_forces(i, j, d, r, u_l, f, do_stress)
910            endif
# Line 833 | Line 912 | contains
912         endif
913      endif
914  
915 <    if (FF_uses_Sticky .and. SimUsesSticky()) then
915 >    if (FF_uses_Sticky .and. SIM_uses_sticky) then
916  
917 <       if ( propertyMapI(3, me_i) .and. propertyMapJ(3, me_j)) then
917 >       if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
918            call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
919                 do_pot, do_stress)
920         endif
921      endif
922  
923  
924 <    if (FF_uses_GB .and. SimUsesGB()) then
924 >    if (FF_uses_GB .and. SIM_uses_GB) then
925        
926 <       if ( propertyMapI(4, me_i) .and. propertyMapJ(4, me_j)) then
926 >       if ( PropertyMap(me_i)%is_GB .and. PropertyMap(me_j)%is_GB) then
927            call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, &
928                 do_pot, do_stress)          
929         endif
# Line 853 | Line 932 | contains
932      
933  
934    
935 <   if (FF_uses_EAM .and. SimUsesEAM()) then
935 >   if (FF_uses_EAM .and. SIM_uses_EAM) then
936        
937 <      if ( propertyMapI(5, me_i) .and. propertyMapJ(5, me_j)) then
937 >      if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
938           call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
939        endif
940  
# Line 900 | Line 979 | contains
979    
980   #endif
981      
982 <   if (FF_uses_EAM .and. SimUsesEAM()) then
983 <      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
984 <      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
906 <      
907 <      if ( is_EAM_i .and. is_EAM_j ) &
982 >   if (FF_uses_EAM .and. SIM_uses_EAM) then
983 >
984 >      if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
985             call calc_EAM_prepair_rho(i, j, d, r, rijsq )
909   endif
986  
987 +   endif
988 +  
989   end subroutine do_prepair
990  
991  
# Line 917 | Line 995 | contains
995      integer :: nlocal
996      real( kind = dp ) :: pot
997  
998 <    if (FF_uses_EAM .and. SimUsesEAM()) then
998 >    if (FF_uses_EAM .and. SIM_uses_EAM) then
999         call calc_EAM_preforce_Frho(nlocal,pot)
1000      endif
1001  
# Line 936 | Line 1014 | contains
1014      d(1:3) = q_j(1:3) - q_i(1:3)
1015  
1016      ! Wrap back into periodic box if necessary
1017 <    if ( SimUsesPBC() ) then
1017 >    if ( SIM_uses_PBC ) then
1018        
1019         if( .not.boxIsOrthorhombic ) then
1020            ! calc the scaled coordinates.
# Line 976 | Line 1054 | contains
1054      
1055    end subroutine get_interatomic_vector
1056    
979  subroutine check_initialization(error)
980    integer, intent(out) :: error
981    
982    error = 0
983    ! Make sure we are properly initialized.
984    if (.not. do_forces_initialized) then
985       write(*,*) "Forces not initialized"
986       error = -1
987       return
988    endif
989
990 #ifdef IS_MPI
991    if (.not. isMPISimSet()) then
992       write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
993       error = -1
994       return
995    endif
996 #endif
997    
998    return
999  end subroutine check_initialization
1000
1001  
1057    subroutine zero_work_arrays()
1058      
1059   #ifdef IS_MPI
# Line 1031 | Line 1086 | contains
1086   #endif
1087  
1088  
1089 <    if (FF_uses_EAM .and. SimUsesEAM()) then
1089 >    if (FF_uses_EAM .and. SIM_uses_EAM) then
1090         call clean_EAM()
1091      endif
1092  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines