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 388 by chuckv, Fri Mar 21 22:11:50 2003 UTC vs.
Revision 841 by mmeineke, Wed Oct 29 17:55:28 2003 UTC

# Line 4 | Line 4
4  
5   !! @author Charles F. Vardeman II
6   !! @author Matthew Meineke
7 < !! @version $Id: do_Forces.F90,v 1.2 2003-03-21 22:11:50 chuckv Exp $, $Date: 2003-03-21 22:11:50 $, $Name: not supported by cvs2svn $, $Revision: 1.2 $
7 > !! @version $Id: do_Forces.F90,v 1.34 2003-10-29 17:55:28 mmeineke Exp $, $Date: 2003-10-29 17:55:28 $, $Name: not supported by cvs2svn $, $Revision: 1.34 $
8  
9   module do_Forces
10    use force_globals
# Line 17 | Line 17 | module do_Forces
17    use dipole_dipole
18    use reaction_field
19    use gb_pair
20 +  use vector_class
21 +  use eam
22 +  use status
23   #ifdef IS_MPI
24    use mpiSimulation
25   #endif
# Line 27 | Line 30 | module do_Forces
30   #define __FORTRAN90
31   #include "fForceField.h"
32  
33 <  logical, save :: do_forces_initialized = .false.
33 >  logical, save :: do_forces_initialized = .false., haveRlist = .false.
34 >  logical, save :: havePolicies = .false.
35    logical, save :: FF_uses_LJ
36    logical, save :: FF_uses_sticky
37    logical, save :: FF_uses_dipoles
# Line 35 | Line 39 | module do_Forces
39    logical, save :: FF_uses_GB
40    logical, save :: FF_uses_EAM
41  
42 +  real(kind=dp), save :: rlist, rlistsq
43 +
44    public :: init_FF
45    public :: do_force_loop
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
54 + #endif
55 +
56   contains
57  
58 +  subroutine setRlistDF( this_rlist )
59 +    
60 +    real(kind=dp) :: this_rlist
61 +
62 +    rlist = this_rlist
63 +    rlistsq = rlist * rlist
64 +    
65 +    haveRlist = .true.
66 +    if( havePolicies ) do_forces_initialized = .true.
67 +
68 +  end subroutine setRlistDF    
69 +
70    subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
71  
72      integer, intent(in) :: LJMIXPOLICY
# Line 87 | Line 114 | contains
114      !! check to make sure the FF_uses_RF setting makes sense
115      
116      if (FF_uses_dipoles) then
90 –       rrf = getRrf()
91 –       rt = getRt()      
92 –       call initialize_dipole(rrf, rt)
117         if (FF_uses_RF) then
118            dielect = getDielect()
119 <          call initialize_rf(rrf, rt, dielect)
119 >          call initialize_rf(dielect)
120         endif
121      else
122         if (FF_uses_RF) then          
# Line 100 | Line 124 | contains
124            thisStat = -1
125            return
126         endif
127 <    endif
127 >    endif
128  
129      if (FF_uses_LJ) then
130        
107 –       call getRcut(rcut)
108 –
131         select case (LJMIXPOLICY)
132         case (LB_MIXING_RULE)
133 <          call init_lj_FF(LB_MIXING_RULE, rcut, my_status)            
133 >          call init_lj_FF(LB_MIXING_RULE, my_status)            
134         case (EXPLICIT_MIXING_RULE)
135 <          call init_lj_FF(EXPLICIT_MIXING_RULE, rcut, my_status)
135 >          call init_lj_FF(EXPLICIT_MIXING_RULE, my_status)
136         case default
137            write(default_error,*) 'unknown LJ Mixing Policy!'
138            thisStat = -1
# Line 129 | Line 151 | contains
151            return
152         end if
153      endif
154 +
155 +
156 +    if (FF_uses_EAM) then
157 +         call init_EAM_FF(my_status)
158 +       if (my_status /= 0) then
159 +          write(*,*) "init_EAM_FF returned a bad status"
160 +          thisStat = -1
161 +          return
162 +       end if
163 +    endif
164 +
165 +
166      
167      if (FF_uses_GB) then
168         call check_gb_pair_FF(my_status)
# Line 140 | Line 174 | contains
174  
175      if (FF_uses_GB .and. FF_uses_LJ) then
176      endif
177 +    if (.not. do_forces_initialized) then
178 +       !! Create neighbor lists
179 +       call expandNeighborList(getNlocal(), my_status)
180 +       if (my_Status /= 0) then
181 +          write(default_error,*) "SimSetup: ExpandNeighborList returned error."
182 +          thisStat = -1
183 +          return
184 +       endif
185 +    endif
186 +    
187  
188 +    havePolicies = .true.
189 +    if( haveRlist ) do_forces_initialized = .true.
190  
145 –    do_forces_initialized = .true.    
146 –
191    end subroutine init_FF
192    
193  
# Line 167 | Line 211 | contains
211      logical ( kind = 2) :: do_pot_c, do_stress_c
212      logical :: do_pot
213      logical :: do_stress
214 < #ifdef IS_MPI
214 > #ifdef IS_MPI
215      real( kind = DP ) :: pot_local
216      integer :: nrow
217      integer :: ncol
218 +    integer :: nprocs
219   #endif
220      integer :: nlocal
221      integer :: natoms    
222      logical :: update_nlist  
223      integer :: i, j, jbeg, jend, jnab
224      integer :: nlist
225 <    real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut
225 >    real( kind = DP ) ::  rijsq
226      real(kind=dp),dimension(3) :: d
227      real(kind=dp) :: rfpot, mu_i, virial
228      integer :: me_i
# Line 186 | Line 231 | contains
231      integer :: listerror, error
232      integer :: localError
233  
234 +    real(kind=dp) :: listSkin = 1.0
235 +    
236 +
237      !! initialize local variables  
238  
239   #ifdef IS_MPI
240 +    pot_local = 0.0_dp
241      nlocal = getNlocal()
242      nrow   = getNrow(plan_row)
243      ncol   = getNcol(plan_col)
# Line 197 | Line 246 | contains
246      natoms = nlocal
247   #endif
248  
200 –    call getRcut(rcut,rc2=rcutsq)
201 –    call getRlist(rlist,rlistsq)
202 –    
249      call check_initialization(localError)
250      if ( localError .ne. 0 ) then
251 +       call handleError("do_force_loop","Not Initialized")
252         error = -1
253         return
254      end if
# Line 210 | Line 257 | contains
257      do_pot = do_pot_c
258      do_stress = do_stress_c
259  
260 +
261      ! Gather all information needed by all force loops:
262      
263   #ifdef IS_MPI    
# Line 226 | Line 274 | contains
274      endif
275      
276   #endif
277 <    
277 >
278 > !! Begin force loop timing:
279 > #ifdef PROFILE
280 >    call cpu_time(forceTimeInitial)
281 >    nloops = nloops + 1
282 > #endif
283 >  
284      if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
285         !! See if we need to update neighbor lists
286 <       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)  
286 >       call checkNeighborList(nlocal, q, listSkin, update_nlist)  
287         !! if_mpi_gather_stuff_for_prepair
288         !! do_prepair_loop_if_needed
289         !! if_mpi_scatter_stuff_from_prepair
290         !! if_mpi_gather_stuff_from_prepair_to_main_loop
291 <    else
292 <       !! See if we need to update neighbor lists
293 <       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)  
291 >    
292 > !--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
293 > #ifdef IS_MPI
294 >    
295 >    if (update_nlist) then
296 >      
297 >       !! save current configuration, construct neighbor list,
298 >       !! and calculate forces
299 >       call saveNeighborList(nlocal, q)
300 >      
301 >       neighborListSize = size(list)
302 >       nlist = 0      
303 >      
304 >       do i = 1, nrow
305 >          point(i) = nlist + 1
306 >          
307 >          prepair_inner: do j = 1, ncol
308 >            
309 >             if (skipThisPair(i,j)) cycle prepair_inner
310 >            
311 >             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
312 >            
313 >             if (rijsq < rlistsq) then            
314 >                
315 >                nlist = nlist + 1
316 >                
317 >                if (nlist > neighborListSize) then
318 >                   call expandNeighborList(nlocal, listerror)
319 >                   if (listerror /= 0) then
320 >                      error = -1
321 >                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
322 >                      return
323 >                   end if
324 >                   neighborListSize = size(list)
325 >                endif
326 >                
327 >                list(nlist) = j
328 >                call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)                      
329 >             endif
330 >          enddo prepair_inner
331 >       enddo
332 >
333 >       point(nrow + 1) = nlist + 1
334 >      
335 >    else  !! (of update_check)
336 >
337 >       ! use the list to find the neighbors
338 >       do i = 1, nrow
339 >          JBEG = POINT(i)
340 >          JEND = POINT(i+1) - 1
341 >          ! check thiat molecule i has neighbors
342 >          if (jbeg .le. jend) then
343 >            
344 >             do jnab = jbeg, jend
345 >                j = list(jnab)
346 >
347 >                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
348 >                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
349 >                     u_l, A, f, t, pot_local)
350 >
351 >             enddo
352 >          endif
353 >       enddo
354      endif
355      
356 + #else
357 +    
358 +    if (update_nlist) then
359 +      
360 +       ! save current configuration, contruct neighbor list,
361 +       ! and calculate forces
362 +       call saveNeighborList(natoms, q)
363 +      
364 +       neighborListSize = size(list)
365 +  
366 +       nlist = 0
367 +
368 +       do i = 1, natoms-1
369 +          point(i) = nlist + 1
370 +          
371 +          prepair_inner: do j = i+1, natoms
372 +            
373 +             if (skipThisPair(i,j))  cycle prepair_inner
374 +                          
375 +             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
376 +          
377 +
378 +             if (rijsq < rlistsq) then
379 +
380 +          
381 +                nlist = nlist + 1
382 +              
383 +                if (nlist > neighborListSize) then
384 +                   call expandNeighborList(natoms, listerror)
385 +                   if (listerror /= 0) then
386 +                      error = -1
387 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
388 +                      return
389 +                   end if
390 +                   neighborListSize = size(list)
391 +                endif
392 +                
393 +                list(nlist) = j
394 +                
395 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
396 +                        u_l, A, f, t, pot)
397 +                
398 +             endif
399 +          enddo prepair_inner
400 +       enddo
401 +      
402 +       point(natoms) = nlist + 1
403 +      
404 +    else !! (update)
405 +  
406 +       ! use the list to find the neighbors
407 +       do i = 1, natoms-1
408 +          JBEG = POINT(i)
409 +          JEND = POINT(i+1) - 1
410 +          ! check thiat molecule i has neighbors
411 +          if (jbeg .le. jend) then
412 +            
413 +             do jnab = jbeg, jend
414 +                j = list(jnab)
415 +
416 +                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
417 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
418 +                     u_l, A, f, t, pot)
419 +
420 +             enddo
421 +          endif
422 +       enddo
423 +    endif    
424 + #endif
425 +    !! Do rest of preforce calculations
426 +    !! do necessary preforce calculations  
427 +    call do_preforce(nlocal,pot)
428 +   ! we have already updated the neighbor list set it to false...
429 +   update_nlist = .false.
430 +    else
431 +       !! See if we need to update neighbor lists for non pre-pair
432 +       call checkNeighborList(nlocal, q, listSkin, update_nlist)  
433 +    endif
434 +
435 +
436 +
437 +
438 +
439 + !---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
440 +
441 +
442 +
443 +
444 +  
445   #ifdef IS_MPI
446      
447      if (update_nlist) then
448        
449         !! save current configuration, construct neighbor list,
450         !! and calculate forces
451 <       call saveNeighborList(q)
451 >       call saveNeighborList(nlocal, q)
452        
453         neighborListSize = size(list)
454         nlist = 0      
# Line 259 | Line 462 | contains
462              
463               call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
464              
465 <             if (rijsq <  rlistsq) then            
465 >             if (rijsq < rlistsq) then            
466                  
467                  nlist = nlist + 1
468                  
# Line 275 | Line 478 | contains
478                  
479                  list(nlist) = j
480                                  
481 <                if (rijsq <  rcutsq) then
482 <                   call do_pair(i, j, rijsq, d, do_pot, do_stress, &
483 <                        u_l, A, f, t,pot)
281 <                endif
481 >                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
482 >                     u_l, A, f, t, pot_local)
483 >                
484               endif
485            enddo inner
486         enddo
# Line 299 | Line 501 | contains
501  
502                  call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
503                  call do_pair(i, j, rijsq, d, do_pot, do_stress, &
504 <                     u_l, A, f, t,pot)
504 >                     u_l, A, f, t, pot_local)
505  
506               enddo
507            endif
# Line 312 | Line 514 | contains
514        
515         ! save current configuration, contruct neighbor list,
516         ! and calculate forces
517 <       call saveNeighborList(q)
517 >       call saveNeighborList(natoms, q)
518        
519         neighborListSize = size(list)
520    
# Line 328 | Line 530 | contains
530               call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
531            
532  
533 <             if (rijsq <  rlistsq) then
533 >             if (rijsq < rlistsq) then
534                  
535                  nlist = nlist + 1
536                
# Line 344 | Line 546 | contains
546                  
547                  list(nlist) = j
548                  
549 <                if (rijsq <  rcutsq) then
550 <                   call do_pair(i, j, rijsq, d, do_pot, do_stress, &
551 <                        u_l, A, f, t,pot)
350 <                endif
549 >                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
550 >                        u_l, A, f, t, pot)
551 >                
552               endif
553            enddo inner
554         enddo
# Line 368 | Line 569 | contains
569  
570                  call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
571                  call do_pair(i, j, rijsq, d, do_pot, do_stress, &
572 <                     u_l, A, f, t,pot)
572 >                     u_l, A, f, t, pot)
573  
574               enddo
575            endif
# Line 378 | Line 579 | contains
579   #endif
580      
581      ! phew, done with main loop.
582 <    
582 >
583 > !! Do timing
584 > #ifdef PROFILE
585 >    call cpu_time(forceTimeFinal)
586 >    forceTime = forceTime + forceTimeFinal - forceTimeInitial
587 > #endif
588 >
589 >
590   #ifdef IS_MPI
591      !!distribute forces
592 <    
593 <    call scatter(f_Row,f,plan_row3d)
592 >  
593 >    f_temp = 0.0_dp
594 >    call scatter(f_Row,f_temp,plan_row3d)
595 >    do i = 1,nlocal
596 >       f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
597 >    end do
598 >
599 >    f_temp = 0.0_dp
600      call scatter(f_Col,f_temp,plan_col3d)
601      do i = 1,nlocal
602         f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
603      end do
604      
605      if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
606 <       call scatter(t_Row,t,plan_row3d)
606 >       t_temp = 0.0_dp
607 >       call scatter(t_Row,t_temp,plan_row3d)
608 >       do i = 1,nlocal
609 >          t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
610 >       end do
611 >       t_temp = 0.0_dp
612         call scatter(t_Col,t_temp,plan_col3d)
613        
614         do i = 1,nlocal
# Line 400 | Line 619 | contains
619      if (do_pot) then
620         ! scatter/gather pot_row into the members of my column
621         call scatter(pot_Row, pot_Temp, plan_row)
622 <      
622 >
623         ! scatter/gather pot_local into all other procs
624         ! add resultant to get total pot
625         do i = 1, nlocal
626            pot_local = pot_local + pot_Temp(i)
627         enddo
628 +      
629 +       pot_Temp = 0.0_DP
630  
410 –       pot_Temp = 0.0_DP
411 –
631         call scatter(pot_Col, pot_Temp, plan_col)
632         do i = 1, nlocal
633            pot_local = pot_local + pot_Temp(i)
634         enddo
635 <      
635 >
636      endif    
637   #endif
638  
# Line 461 | Line 680 | contains
680   #ifdef IS_MPI
681  
682      if (do_pot) then
683 <       pot = pot_local
683 >       pot = pot + pot_local
684         !! we assume the c code will do the allreduce to get the total potential
685         !! we could do it right here if we needed to...
686      endif
687  
688      if (do_stress) then
689 <       call mpi_allreduce(tau, tau_Temp,9,mpi_double_precision,mpi_sum, &
689 >      call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
690              mpi_comm_world,mpi_err)
691 <       call mpi_allreduce(virial, virial_Temp,1,mpi_double_precision,mpi_sum, &
691 >       call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
692              mpi_comm_world,mpi_err)
693      endif
694  
# Line 481 | Line 700 | contains
700      endif
701  
702   #endif
703 <    
703 >
704 > #ifdef PROFILE
705 >    if (do_pot) then
706 >
707 > #ifdef IS_MPI
708 >
709 >      
710 >       call printCommTime()
711 >
712 >       call mpi_allreduce(forceTime,globalForceTime,1,MPI_DOUBLE_PRECISION, &
713 >            mpi_sum,mpi_comm_world,mpi_err)
714 >
715 >       call mpi_allreduce(forceTime,maxForceTime,1,MPI_DOUBLE_PRECISION, &
716 >            MPI_MAX,mpi_comm_world,mpi_err)
717 >      
718 >       call mpi_comm_size( MPI_COMM_WORLD, nprocs,mpi_err)
719 >      
720 >       if (getMyNode() == 0) then
721 >          write(*,*) "Total processor time spent in force calculations is: ", globalForceTime
722 >          write(*,*) "Total Time spent in force loop per processor is: ", globalforceTime/nprocs
723 >          write(*,*) "Maximum force time on any processor is: ", maxForceTime
724 >       end if
725 > #else
726 >       write(*,*) "Time spent in force loop is: ", forceTime
727 > #endif
728 >
729 >    
730 >    endif
731 >
732 > #endif
733 >
734 >
735 >    write(*,*) "In F, after: Tau = "
736 >    write(*,*) tau(1), tau(2), tau(3)
737 >    write(*,*) tau(4), tau(5), tau(6)
738 >    write(*,*) tau(7), tau(8), tau(9)
739 >    write(*,*) " "
740 >
741 >
742    end subroutine do_force_loop
743  
744 <  subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t,pot)
744 >  subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
745  
746      real( kind = dp ) :: pot
747 <    real( kind = dp ), dimension(:,:) :: u_l
748 <    real (kind=dp), dimension(:,:) :: A
749 <    real (kind=dp), dimension(:,:) :: f
750 <    real (kind=dp), dimension(:,:) :: t
747 >    real( kind = dp ), dimension(3,getNlocal()) :: u_l
748 >    real (kind=dp), dimension(9,getNlocal()) :: A
749 >    real (kind=dp), dimension(3,getNlocal()) :: f
750 >    real (kind=dp), dimension(3,getNlocal()) :: t
751  
752      logical, intent(inout) :: do_pot, do_stress
753      integer, intent(in) :: i, j
# Line 500 | Line 757 | contains
757      logical :: is_LJ_i, is_LJ_j
758      logical :: is_DP_i, is_DP_j
759      logical :: is_GB_i, is_GB_j
760 +    logical :: is_EAM_i,is_EAM_j
761      logical :: is_Sticky_i, is_Sticky_j
762      integer :: me_i, me_j
763  
764      r = sqrt(rijsq)
765  
766   #ifdef IS_MPI
767 +    if (tagRow(i) .eq. tagColumn(j)) then
768 +       write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
769 +    endif
770  
771      me_i = atid_row(i)
772      me_j = atid_col(j)
# Line 530 | Line 791 | contains
791         call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
792        
793         if ( is_DP_i .and. is_DP_j ) then
533 –          
794            call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
795                 do_pot, do_stress)
796            if (FF_uses_RF .and. SimUsesRF()) then
# Line 555 | Line 815 | contains
815  
816      if (FF_uses_GB .and. SimUsesGB()) then
817  
818 +
819         call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
820         call getElementProperty(atypes, me_j, "is_GB", is_GB_j)
821        
# Line 564 | Line 825 | contains
825         endif
826      endif
827      
828 +
829 +  
830 +   if (FF_uses_EAM .and. SimUsesEAM()) then
831 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
832 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
833 +      
834 +      if ( is_EAM_i .and. is_EAM_j ) &
835 +           call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
836 +   endif
837 +
838 +
839 +
840 +
841    end subroutine do_pair
842  
843  
844 +
845 +  subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
846 +   real( kind = dp ) :: pot
847 +   real( kind = dp ), dimension(3,getNlocal()) :: u_l
848 +   real (kind=dp), dimension(9,getNlocal()) :: A
849 +   real (kind=dp), dimension(3,getNlocal()) :: f
850 +   real (kind=dp), dimension(3,getNlocal()) :: t
851 +  
852 +   logical, intent(inout) :: do_pot, do_stress
853 +   integer, intent(in) :: i, j
854 +   real ( kind = dp ), intent(inout)    :: rijsq
855 +   real ( kind = dp )                :: r
856 +   real ( kind = dp ), intent(inout) :: d(3)
857 +  
858 +   logical :: is_EAM_i, is_EAM_j
859 +  
860 +   integer :: me_i, me_j
861 +  
862 +   r = sqrt(rijsq)
863 +  
864 +
865 + #ifdef IS_MPI
866 +   if (tagRow(i) .eq. tagColumn(j)) then
867 +      write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
868 +   endif
869 +  
870 +   me_i = atid_row(i)
871 +   me_j = atid_col(j)
872 +  
873 + #else
874 +  
875 +   me_i = atid(i)
876 +   me_j = atid(j)
877 +  
878 + #endif
879 +    
880 +   if (FF_uses_EAM .and. SimUsesEAM()) then
881 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
882 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
883 +      
884 +      if ( is_EAM_i .and. is_EAM_j ) &
885 +           call calc_EAM_prepair_rho(i, j, d, r, rijsq )
886 +   endif
887 +
888 + end subroutine do_prepair
889 +
890 +
891 +
892 +
893 +  subroutine do_preforce(nlocal,pot)
894 +    integer :: nlocal
895 +    real( kind = dp ) :: pot
896 +
897 +    if (FF_uses_EAM .and. SimUsesEAM()) then
898 +       call calc_EAM_preforce_Frho(nlocal,pot)
899 +    endif
900 +
901 +
902 +  end subroutine do_preforce
903 +  
904 +  
905    subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
906      
907      real (kind = dp), dimension(3) :: q_i
908      real (kind = dp), dimension(3) :: q_j
909      real ( kind = dp ), intent(out) :: r_sq
910 <    real( kind = dp ) :: d(3)
911 <    real( kind = dp ) :: d_old(3)
912 <    d(1:3) = q_i(1:3) - q_j(1:3)
913 <    d_old = d
910 >    real( kind = dp ) :: d(3), scaled(3)
911 >    integer i
912 >
913 >    d(1:3) = q_j(1:3) - q_i(1:3)
914 >
915      ! Wrap back into periodic box if necessary
916      if ( SimUsesPBC() ) then
917 +      
918 +       if( .not.boxIsOrthorhombic ) then
919 +          ! calc the scaled coordinates.
920 +          
921 +          scaled = matmul(HmatInv, d)
922 +          
923 +          ! wrap the scaled coordinates
924  
925 <       d(1:3) = d(1:3) - box(1:3) * sign(1.0_dp,d(1:3)) * &
926 <            int(abs(d(1:3)/box(1:3)) + 0.5_dp)
925 >          scaled = scaled  - anint(scaled)
926 >          
927  
928 +          ! calc the wrapped real coordinates from the wrapped scaled
929 +          ! coordinates
930 +
931 +          d = matmul(Hmat,scaled)
932 +
933 +       else
934 +          ! calc the scaled coordinates.
935 +          
936 +          do i = 1, 3
937 +             scaled(i) = d(i) * HmatInv(i,i)
938 +            
939 +             ! wrap the scaled coordinates
940 +            
941 +             scaled(i) = scaled(i) - anint(scaled(i))
942 +            
943 +             ! calc the wrapped real coordinates from the wrapped scaled
944 +             ! coordinates
945 +
946 +             d(i) = scaled(i)*Hmat(i,i)
947 +          enddo
948 +       endif
949 +      
950      endif
951 +    
952      r_sq = dot_product(d,d)
953 <        
953 >    
954    end subroutine get_interatomic_vector
955 <
955 >  
956    subroutine check_initialization(error)
957      integer, intent(out) :: error
958      
959      error = 0
960      ! Make sure we are properly initialized.
961      if (.not. do_forces_initialized) then
962 +       write(*,*) "Forces not initialized"
963         error = -1
964         return
965      endif
# Line 640 | Line 1007 | contains
1007  
1008   #endif
1009  
1010 +
1011 +    if (FF_uses_EAM .and. SimUsesEAM()) then
1012 +       call clean_EAM()
1013 +    endif
1014 +
1015 +
1016 +
1017 +
1018 +
1019      rf = 0.0_dp
1020      tau_Temp = 0.0_dp
1021      virial_Temp = 0.0_dp
646 –    
1022    end subroutine zero_work_arrays
1023    
1024    function skipThisPair(atom1, atom2) result(skip_it)
# Line 687 | Line 1062 | contains
1062   #else
1063      unique_id_2 = atom2
1064   #endif
1065 <    
1065 >
1066   #ifdef IS_MPI
1067      !! this situation should only arise in MPI simulations
1068      if (unique_id_1 == unique_id_2) then
# Line 697 | Line 1072 | contains
1072      
1073      !! this prevents us from doing the pair on multiple processors
1074      if (unique_id_1 < unique_id_2) then
1075 <       if (mod(unique_id_1 + unique_id_2,2) == 0) skip_it = .true.
1076 <       return
1075 >       if (mod(unique_id_1 + unique_id_2,2) == 0) then
1076 >          skip_it = .true.
1077 >          return
1078 >       endif
1079      else                
1080 <       if (mod(unique_id_1 + unique_id_2,2) == 1) skip_it = .true.
1081 <       return
1080 >       if (mod(unique_id_1 + unique_id_2,2) == 1) then
1081 >          skip_it = .true.
1082 >          return
1083 >       endif
1084      endif
1085   #endif
1086 <
1086 >
1087      !! the rest of these situations can happen in all simulations:
1088      do i = 1, nExcludes_global      
1089         if ((excludesGlobal(i) == unique_id_1) .or. &
# Line 713 | Line 1092 | contains
1092            return
1093         endif
1094      enddo
1095 <
1095 >
1096      do i = 1, nExcludes_local
1097         if (excludesLocal(1,i) == unique_id_1) then
1098            if (excludesLocal(2,i) == unique_id_2) then
# Line 749 | Line 1128 | end module do_Forces
1128      doesit = FF_uses_RF
1129    end function FF_RequiresPostpairCalc
1130    
1131 + !! This cleans componets of force arrays belonging only to fortran
1132 +
1133   end module do_Forces

Diff Legend

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