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 726 by tim, Tue Aug 26 20:37:30 2003 UTC vs.
Revision 898 by chuckv, Mon Jan 5 22:49:14 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.30 2003-08-26 20:37:30 tim Exp $, $Date: 2003-08-26 20:37:30 $, $Name: not supported by cvs2svn $, $Revision: 1.30 $
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 $
8  
9   module do_Forces
10    use force_globals
# Line 46 | Line 46 | module do_Forces
46    public :: setRlistDF
47  
48   #ifdef PROFILE
49 <  real(kind = dp) :: forceTime
50 <  real(kind = dp) :: forceTimeInitial, forceTimeFinal
51 <  real(kind = dp) :: globalForceTime
52 <  real(kind = dp) :: maxForceTime
53 <  integer, save :: nloops = 0
49 >  public :: getforcetime
50 >  real, save :: forceTime = 0
51 >  real :: forceTimeInitial, forceTimeFinal
52 >  integer :: nLoops
53   #endif
54  
55 +  logical, allocatable :: propertyMapI(:,:)
56 +  logical, allocatable :: propertyMapJ(:,:)
57 +
58   contains
59  
60    subroutine setRlistDF( this_rlist )
# Line 154 | Line 156 | contains
156  
157  
158      if (FF_uses_EAM) then
159 <       call init_EAM_FF(my_status)
159 >         call init_EAM_FF(my_status)
160         if (my_status /= 0) then
161 +          write(*,*) "init_EAM_FF returned a bad status"
162            thisStat = -1
163            return
164         end if
# Line 175 | Line 178 | contains
178      endif
179      if (.not. do_forces_initialized) then
180         !! Create neighbor lists
181 <       call expandNeighborList(getNlocal(), my_status)
181 >       call expandNeighborList(nLocal, my_status)
182         if (my_Status /= 0) then
183            write(default_error,*) "SimSetup: ExpandNeighborList returned error."
184            thisStat = -1
# Line 195 | Line 198 | contains
198    subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, &
199         error)
200      !! Position array provided by C, dimensioned by getNlocal
201 <    real ( kind = dp ), dimension(3,getNlocal()) :: q
201 >    real ( kind = dp ), dimension(3,nLocal) :: q
202      !! Rotation Matrix for each long range particle in simulation.
203 <    real( kind = dp), dimension(9,getNlocal()) :: A    
203 >    real( kind = dp), dimension(9,nLocal) :: A    
204      !! Unit vectors for dipoles (lab frame)
205 <    real( kind = dp ), dimension(3,getNlocal()) :: u_l
205 >    real( kind = dp ), dimension(3,nLocal) :: u_l
206      !! Force array provided by C, dimensioned by getNlocal
207 <    real ( kind = dp ), dimension(3,getNlocal()) :: f
207 >    real ( kind = dp ), dimension(3,nLocal) :: f
208      !! Torsion array provided by C, dimensioned by getNlocal
209 <    real( kind = dp ), dimension(3,getNlocal()) :: t    
209 >    real( kind = dp ), dimension(3,nLocal) :: t    
210 >
211      !! Stress Tensor
212      real( kind = dp), dimension(9) :: tau  
213      real ( kind = dp ) :: pot
# Line 216 | Line 220 | contains
220      integer :: ncol
221      integer :: nprocs
222   #endif
219 –    integer :: nlocal
223      integer :: natoms    
224      logical :: update_nlist  
225      integer :: i, j, jbeg, jend, jnab
# Line 224 | Line 227 | contains
227      real( kind = DP ) ::  rijsq
228      real(kind=dp),dimension(3) :: d
229      real(kind=dp) :: rfpot, mu_i, virial
230 <    integer :: me_i
230 >    integer :: me_i, me_j
231      logical :: is_dp_i
232      integer :: neighborListSize
233      integer :: listerror, error
234      integer :: localError
235 +    integer :: propPack_i, propPack_j
236  
237 <    real(kind=dp) :: listSkin = 1.0
234 <    
237 >    real(kind=dp) :: listSkin = 1.0  
238  
239      !! initialize local variables  
240  
241   #ifdef IS_MPI
242      pot_local = 0.0_dp
240 –    nlocal = getNlocal()
243      nrow   = getNrow(plan_row)
244      ncol   = getNcol(plan_col)
245   #else
244 –    nlocal = getNlocal()
246      natoms = nlocal
247   #endif
248  
# Line 256 | Line 257 | contains
257      do_pot = do_pot_c
258      do_stress = do_stress_c
259  
260 <
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
300 >
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 >
331      ! Gather all information needed by all force loops:
332      
333   #ifdef IS_MPI    
# Line 444 | Line 515 | contains
515   #ifdef IS_MPI
516      
517      if (update_nlist) then
447 –      
518         !! save current configuration, construct neighbor list,
519         !! and calculate forces
520         call saveNeighborList(nlocal, q)
# Line 453 | Line 523 | contains
523         nlist = 0      
524        
525         do i = 1, nrow
526 +
527            point(i) = nlist + 1
528            
529            inner: do j = 1, ncol
# Line 510 | Line 581 | contains
581   #else
582      
583      if (update_nlist) then
584 <      
584 >
585         ! save current configuration, contruct neighbor list,
586         ! and calculate forces
587         call saveNeighborList(natoms, q)
# Line 647 | Line 718 | contains
718            end do
719   #endif
720            
721 <          do i = 1, getNlocal()
721 >          do i = 1, nLocal
722  
723               rfpot = 0.0_DP
724   #ifdef IS_MPI
# Line 697 | Line 768 | contains
768         tau = tau_Temp
769         virial = virial_Temp
770      endif
771 <
771 >    
772   #endif
773 <
774 < #ifdef PROFILE
775 <    if (do_pot) then
705 <
706 < #ifdef IS_MPI
707 <
708 <      
709 <       call printCommTime()
710 <
711 <       call mpi_allreduce(forceTime,globalForceTime,1,MPI_DOUBLE_PRECISION, &
712 <            mpi_sum,mpi_comm_world,mpi_err)
713 <
714 <       call mpi_allreduce(forceTime,maxForceTime,1,MPI_DOUBLE_PRECISION, &
715 <            MPI_MAX,mpi_comm_world,mpi_err)
716 <      
717 <       call mpi_comm_size( MPI_COMM_WORLD, nprocs,mpi_err)
718 <      
719 <       if (getMyNode() == 0) then
720 <          write(*,*) "Total processor time spent in force calculations is: ", globalForceTime
721 <          write(*,*) "Total Time spent in force loop per processor is: ", globalforceTime/nprocs
722 <          write(*,*) "Maximum force time on any processor is: ", maxForceTime
723 <       end if
724 < #else
725 <       write(*,*) "Time spent in force loop is: ", forceTime
726 < #endif
727 <
728 <    
729 <    endif
730 <
731 < #endif
732 <
733 <
734 <
773 >    
774 >    
775 >    
776    end subroutine do_force_loop
777  
778    subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
779  
780      real( kind = dp ) :: pot
781 <    real( kind = dp ), dimension(3,getNlocal()) :: u_l
782 <    real (kind=dp), dimension(9,getNlocal()) :: A
783 <    real (kind=dp), dimension(3,getNlocal()) :: f
784 <    real (kind=dp), dimension(3,getNlocal()) :: t
781 >    real( kind = dp ), dimension(3,nLocal) :: u_l
782 >    real (kind=dp), dimension(9,nLocal) :: A
783 >    real (kind=dp), dimension(3,nLocal) :: f
784 >    real (kind=dp), dimension(3,nLocal) :: t
785  
786      logical, intent(inout) :: do_pot, do_stress
787      integer, intent(in) :: i, j
# Line 753 | Line 794 | contains
794      logical :: is_EAM_i,is_EAM_j
795      logical :: is_Sticky_i, is_Sticky_j
796      integer :: me_i, me_j
797 <
797 >    integer :: propPack_i
798 >    integer :: propPack_j
799      r = sqrt(rijsq)
800  
801   #ifdef IS_MPI
# Line 770 | Line 812 | contains
812      me_j = atid(j)
813  
814   #endif
815 <
815 >    
816      if (FF_uses_LJ .and. SimUsesLJ()) then
775 –       call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
776 –       call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
817  
818 <       if ( is_LJ_i .and. is_LJ_j ) &
818 >       if ( propertyMapI(1, me_i) .and. propertyMapJ(1, me_j) ) &
819              call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
820 +
821      endif
822  
823      if (FF_uses_dipoles .and. SimUsesDipoles()) then
783 –       call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
784 –       call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
824        
825 <       if ( is_DP_i .and. is_DP_j ) then
825 >       if ( propertyMapI(2, me_i) .and. propertyMapJ(2, me_j)) then
826            call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
827                 do_pot, do_stress)
828            if (FF_uses_RF .and. SimUsesRF()) then
# Line 796 | Line 835 | contains
835  
836      if (FF_uses_Sticky .and. SimUsesSticky()) then
837  
838 <       call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
800 <       call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
801 <
802 <       if ( is_Sticky_i .and. is_Sticky_j ) then
838 >       if ( propertyMapI(3, me_i) .and. propertyMapJ(3, me_j)) then
839            call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
840                 do_pot, do_stress)
841         endif
# Line 807 | Line 843 | contains
843  
844  
845      if (FF_uses_GB .and. SimUsesGB()) then
810 –
811 –
812 –       call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
813 –       call getElementProperty(atypes, me_j, "is_GB", is_GB_j)
846        
847 <       if ( is_GB_i .and. is_GB_j ) then
847 >       if ( propertyMapI(4, me_i) .and. propertyMapJ(4, me_j)) then
848            call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, &
849                 do_pot, do_stress)          
850         endif
851 +
852      endif
853      
854  
855    
856     if (FF_uses_EAM .and. SimUsesEAM()) then
824 –      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
825 –      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
857        
858 <      if ( is_EAM_i .and. is_EAM_j ) &
859 <           call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
858 >      if ( propertyMapI(5, me_i) .and. propertyMapJ(5, me_j)) then
859 >         call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
860 >      endif
861 >
862     endif
830 –
863  
832 –
833 –
864    end subroutine do_pair
865  
866  
867  
868    subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
869     real( kind = dp ) :: pot
870 <   real( kind = dp ), dimension(3,getNlocal()) :: u_l
871 <   real (kind=dp), dimension(9,getNlocal()) :: A
872 <   real (kind=dp), dimension(3,getNlocal()) :: f
873 <   real (kind=dp), dimension(3,getNlocal()) :: t
870 >   real( kind = dp ), dimension(3,nLocal) :: u_l
871 >   real (kind=dp), dimension(9,nLocal) :: A
872 >   real (kind=dp), dimension(3,nLocal) :: f
873 >   real (kind=dp), dimension(3,nLocal) :: t
874    
875     logical, intent(inout) :: do_pot, do_stress
876     integer, intent(in) :: i, j
# Line 1121 | Line 1151 | contains
1151      doesit = FF_uses_RF
1152    end function FF_RequiresPostpairCalc
1153    
1154 + #ifdef PROFILE
1155 +  function getforcetime() result(totalforcetime)
1156 +    real(kind=dp) :: totalforcetime
1157 +    totalforcetime = forcetime
1158 +  end function getforcetime
1159 + #endif
1160 +
1161   !! This cleans componets of force arrays belonging only to fortran
1162  
1163   end module do_Forces

Diff Legend

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