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 462 by gezelter, Sat Apr 5 02:56:27 2003 UTC vs.
Revision 657 by chuckv, Wed Jul 30 21:17:01 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.10 2003-04-05 02:56:27 gezelter Exp $, $Date: 2003-04-05 02:56:27 $, $Name: not supported by cvs2svn $, $Revision: 1.10 $
7 > !! @version $Id: do_Forces.F90,v 1.26 2003-07-30 21:17:01 chuckv Exp $, $Date: 2003-07-30 21:17:01 $, $Name: not supported by cvs2svn $, $Revision: 1.26 $
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   contains
49  
50 +  subroutine setRlistDF( this_rlist )
51 +    
52 +    real(kind=dp) :: this_rlist
53 +
54 +    rlist = this_rlist
55 +    rlistsq = rlist * rlist
56 +    
57 +    haveRlist = .true.
58 +    if( havePolicies ) do_forces_initialized = .true.
59 +
60 +  end subroutine setRlistDF    
61 +
62    subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
63  
64      integer, intent(in) :: LJMIXPOLICY
# Line 87 | Line 106 | contains
106      !! check to make sure the FF_uses_RF setting makes sense
107      
108      if (FF_uses_dipoles) then
90 –       rrf = getRrf()
91 –       rt = getRt()      
92 –       call initialize_dipole(rrf, rt)
109         if (FF_uses_RF) then
110            dielect = getDielect()
111 <          call initialize_rf(rrf, rt, dielect)
111 >          call initialize_rf(dielect)
112         endif
113      else
114         if (FF_uses_RF) then          
# Line 100 | Line 116 | contains
116            thisStat = -1
117            return
118         endif
119 <    endif
119 >    endif
120  
121      if (FF_uses_LJ) then
122        
107 –       call getRcut(rcut)
108 –
123         select case (LJMIXPOLICY)
124         case (LB_MIXING_RULE)
125 <          call init_lj_FF(LB_MIXING_RULE, rcut, my_status)            
125 >          call init_lj_FF(LB_MIXING_RULE, my_status)            
126         case (EXPLICIT_MIXING_RULE)
127 <          call init_lj_FF(EXPLICIT_MIXING_RULE, rcut, my_status)
127 >          call init_lj_FF(EXPLICIT_MIXING_RULE, my_status)
128         case default
129            write(default_error,*) 'unknown LJ Mixing Policy!'
130            thisStat = -1
# Line 129 | Line 143 | contains
143            return
144         end if
145      endif
146 +
147 +
148 +    if (FF_uses_EAM) then
149 +       call init_EAM_FF(my_status)
150 +       if (my_status /= 0) then
151 +          thisStat = -1
152 +          return
153 +       end if
154 +    endif
155 +
156 +
157      
158      if (FF_uses_GB) then
159         call check_gb_pair_FF(my_status)
# Line 140 | Line 165 | contains
165  
166      if (FF_uses_GB .and. FF_uses_LJ) then
167      endif
168 +    if (.not. do_forces_initialized) then
169 +       !! Create neighbor lists
170 +       call expandNeighborList(getNlocal(), my_status)
171 +       if (my_Status /= 0) then
172 +          write(default_error,*) "SimSetup: ExpandNeighborList returned error."
173 +          thisStat = -1
174 +          return
175 +       endif
176 +    endif
177 +    
178  
179 +    havePolicies = .true.
180 +    if( haveRlist ) do_forces_initialized = .true.
181  
145 –    do_forces_initialized = .true.    
146 –
182    end subroutine init_FF
183    
184  
# Line 177 | Line 212 | contains
212      logical :: update_nlist  
213      integer :: i, j, jbeg, jend, jnab
214      integer :: nlist
215 <    real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut
215 >    real( kind = DP ) ::  rijsq
216      real(kind=dp),dimension(3) :: d
217      real(kind=dp) :: rfpot, mu_i, virial
218      integer :: me_i
# Line 186 | Line 221 | contains
221      integer :: listerror, error
222      integer :: localError
223  
224 +    real(kind=dp) :: listSkin = 1.0
225 +    
226 +
227      !! initialize local variables  
228  
229   #ifdef IS_MPI
# Line 197 | Line 235 | contains
235      nlocal = getNlocal()
236      natoms = nlocal
237   #endif
238 <  
201 <    call getRcut(rcut,rc2=rcutsq)
202 <    call getRlist(rlist,rlistsq)
203 <    
238 >    write(*,*) "Inside do_Force Loop"
239      call check_initialization(localError)
240      if ( localError .ne. 0 ) then
241 +       call handleError("do_force_loop","Not Initialized")
242         error = -1
243         return
244      end if
# Line 230 | Line 266 | contains
266      
267      if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
268         !! See if we need to update neighbor lists
269 <       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)  
269 >       call checkNeighborList(nlocal, q, listSkin, update_nlist)  
270         !! if_mpi_gather_stuff_for_prepair
271         !! do_prepair_loop_if_needed
272         !! if_mpi_scatter_stuff_from_prepair
273         !! if_mpi_gather_stuff_from_prepair_to_main_loop
274 <    else
275 <       !! See if we need to update neighbor lists
276 <       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)  
274 >
275 > !--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
276 > #ifdef IS_MPI
277 >    
278 >    if (update_nlist) then
279 >      
280 >       !! save current configuration, construct neighbor list,
281 >       !! and calculate forces
282 >       call saveNeighborList(nlocal, q)
283 >      
284 >       neighborListSize = size(list)
285 >       nlist = 0      
286 >      
287 >       do i = 1, nrow
288 >          point(i) = nlist + 1
289 >          
290 >          prepair_inner: do j = 1, ncol
291 >            
292 >             if (skipThisPair(i,j)) cycle prepair_inner
293 >            
294 >             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
295 >            
296 >             if (rijsq < rlistsq) then            
297 >                
298 >                nlist = nlist + 1
299 >                
300 >                if (nlist > neighborListSize) then
301 >                   call expandNeighborList(nlocal, listerror)
302 >                   if (listerror /= 0) then
303 >                      error = -1
304 >                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
305 >                      return
306 >                   end if
307 >                   neighborListSize = size(list)
308 >                endif
309 >                
310 >                list(nlist) = j
311 >                call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)                      
312 >             endif
313 >          enddo prepair_inner
314 >       enddo
315 >
316 >       point(nrow + 1) = nlist + 1
317 >      
318 >    else  !! (of update_check)
319 >
320 >       ! use the list to find the neighbors
321 >       do i = 1, nrow
322 >          JBEG = POINT(i)
323 >          JEND = POINT(i+1) - 1
324 >          ! check thiat molecule i has neighbors
325 >          if (jbeg .le. jend) then
326 >            
327 >             do jnab = jbeg, jend
328 >                j = list(jnab)
329 >
330 >                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
331 >                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
332 >                     u_l, A, f, t, pot_local)
333 >
334 >             enddo
335 >          endif
336 >       enddo
337      endif
338      
339 + #else
340 +    
341 +    if (update_nlist) then
342 +      
343 +       ! save current configuration, contruct neighbor list,
344 +       ! and calculate forces
345 +       call saveNeighborList(natoms, q)
346 +      
347 +       neighborListSize = size(list)
348 +  
349 +       nlist = 0
350 +      
351 +       do i = 1, natoms-1
352 +          point(i) = nlist + 1
353 +          
354 +          prepair_inner: do j = i+1, natoms
355 +            
356 +             if (skipThisPair(i,j))  cycle prepair_inner
357 +                          
358 +             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
359 +          
360 +
361 +             if (rijsq < rlistsq) then
362 +                
363 +                nlist = nlist + 1
364 +              
365 +                if (nlist > neighborListSize) then
366 +                   call expandNeighborList(natoms, listerror)
367 +                   if (listerror /= 0) then
368 +                      error = -1
369 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
370 +                      return
371 +                   end if
372 +                   neighborListSize = size(list)
373 +                endif
374 +                
375 +                list(nlist) = j
376 +                
377 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
378 +                        u_l, A, f, t, pot)
379 +                
380 +             endif
381 +          enddo prepair_inner
382 +       enddo
383 +      
384 +       point(natoms) = nlist + 1
385 +      
386 +    else !! (update)
387 +      
388 +       ! use the list to find the neighbors
389 +       do i = 1, natoms-1
390 +          JBEG = POINT(i)
391 +          JEND = POINT(i+1) - 1
392 +          ! check thiat molecule i has neighbors
393 +          if (jbeg .le. jend) then
394 +            
395 +             do jnab = jbeg, jend
396 +                j = list(jnab)
397 +
398 +                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
399 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
400 +                     u_l, A, f, t, pot)
401 +
402 +             enddo
403 +          endif
404 +       enddo
405 +    endif    
406 + #endif
407 +    !! Do rest of preforce calculations
408 +   call do_preforce(nlocal,pot)
409 +    else
410 +       !! See if we need to update neighbor lists for non pre-pair
411 +       call checkNeighborList(nlocal, q, listSkin, update_nlist)  
412 +    endif
413 +
414 +
415 +
416 +
417 +
418 + !---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
419 +
420 +
421 +
422 +
423 +  
424   #ifdef IS_MPI
425      
426      if (update_nlist) then
# Line 260 | Line 441 | contains
441              
442               call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
443              
444 <             if (rijsq <  rlistsq) then            
444 >             if (rijsq < rlistsq) then            
445                  
446                  nlist = nlist + 1
447                  
# Line 276 | Line 457 | contains
457                  
458                  list(nlist) = j
459                                  
460 <                if (rijsq <  rcutsq) then
461 <                   call do_pair(i, j, rijsq, d, do_pot, do_stress, &
462 <                        u_l, A, f, t, pot_local)
282 <                endif
460 >                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
461 >                     u_l, A, f, t, pot_local)
462 >                
463               endif
464            enddo inner
465         enddo
# Line 329 | Line 509 | contains
509               call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
510            
511  
512 <             if (rijsq <  rlistsq) then
512 >             if (rijsq < rlistsq) then
513                  
514                  nlist = nlist + 1
515                
# Line 345 | Line 525 | contains
525                  
526                  list(nlist) = j
527                  
528 <                if (rijsq <  rcutsq) then
349 <                   call do_pair(i, j, rijsq, d, do_pot, do_stress, &
528 >                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
529                          u_l, A, f, t, pot)
530 <                endif
530 >                
531               endif
532            enddo inner
533         enddo
# Line 479 | Line 658 | contains
658      endif
659  
660      if (do_stress) then
661 <       call mpi_allreduce(tau, tau_Temp,9,mpi_double_precision,mpi_sum, &
661 >      call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
662              mpi_comm_world,mpi_err)
663 <       call mpi_allreduce(virial, virial_Temp,1,mpi_double_precision,mpi_sum, &
663 >       call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
664              mpi_comm_world,mpi_err)
665      endif
666  
# Line 512 | Line 691 | contains
691      logical :: is_LJ_i, is_LJ_j
692      logical :: is_DP_i, is_DP_j
693      logical :: is_GB_i, is_GB_j
694 +    logical :: is_EAM_i,is_EAM_j
695      logical :: is_Sticky_i, is_Sticky_j
696      integer :: me_i, me_j
697  
698      r = sqrt(rijsq)
699  
700   #ifdef IS_MPI
701 +    if (tagRow(i) .eq. tagColumn(j)) then
702 +       write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
703 +    endif
704  
705      me_i = atid_row(i)
706      me_j = atid_col(j)
# Line 576 | Line 759 | contains
759         endif
760      endif
761      
762 +
763 +  
764 +   if (FF_uses_EAM .and. SimUsesEAM()) then
765 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
766 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
767 +      
768 +      if ( is_EAM_i .and. is_EAM_j ) &
769 +           call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
770 +   endif
771 +
772 +
773 +
774 +
775    end subroutine do_pair
776  
777  
778 +
779 +  subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
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
785 +  
786 +   logical, intent(inout) :: do_pot, do_stress
787 +   integer, intent(in) :: i, j
788 +   real ( kind = dp ), intent(inout)    :: rijsq
789 +   real ( kind = dp )                :: r
790 +   real ( kind = dp ), intent(inout) :: d(3)
791 +  
792 +   logical :: is_EAM_i, is_EAM_j
793 +  
794 +   integer :: me_i, me_j
795 +  
796 +   r = sqrt(rijsq)
797 +  
798 + #ifdef IS_MPI
799 +   if (tagRow(i) .eq. tagColumn(j)) then
800 +      write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
801 +   endif
802 +  
803 +   me_i = atid_row(i)
804 +   me_j = atid_col(j)
805 +  
806 + #else
807 +  
808 +   me_i = atid(i)
809 +   me_j = atid(j)
810 +  
811 + #endif
812 +  
813 +   if (FF_uses_EAM .and. SimUsesEAM()) then
814 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
815 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
816 +      
817 +      if ( is_EAM_i .and. is_EAM_j ) &
818 +           call calc_EAM_prepair_rho(i, j, d, r, rijsq )
819 +   endif
820 +  end subroutine do_prepair
821 +
822 +
823 +
824 +
825 +  subroutine do_preforce(nlocal,pot)
826 +    integer :: nlocal
827 +    real( kind = dp ) :: pot
828 +
829 +   if (FF_uses_EAM .and. SimUsesEAM()) then
830 +      call calc_EAM_preforce_Frho(nlocal,pot)
831 +   endif
832 +
833 +
834 +  end subroutine do_preforce
835 +  
836 +  
837    subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
838      
839      real (kind = dp), dimension(3) :: q_i
840      real (kind = dp), dimension(3) :: q_j
841      real ( kind = dp ), intent(out) :: r_sq
842 <    real( kind = dp ) :: d(3)
843 <    real( kind = dp ) :: d_old(3)
844 <    d(1:3) = q_i(1:3) - q_j(1:3)
845 <    d_old = d
842 >    real( kind = dp ) :: d(3), scaled(3)
843 >    integer i
844 >
845 >    d(1:3) = q_j(1:3) - q_i(1:3)
846 >
847      ! Wrap back into periodic box if necessary
848      if ( SimUsesPBC() ) then
849        
850 <       d(1:3) = d(1:3) - box(1:3) * sign(1.0_dp,d(1:3)) * &
851 <            int(abs(d(1:3)/box(1:3)) + 0.5_dp)
850 >       if( .not.boxIsOrthorhombic ) then
851 >          ! calc the scaled coordinates.
852 >          
853 >          scaled = matmul(HmatInv, d)
854 >          
855 >          ! wrap the scaled coordinates
856 >
857 >          scaled = scaled  - anint(scaled)
858 >          
859 >
860 >          ! calc the wrapped real coordinates from the wrapped scaled
861 >          ! coordinates
862 >
863 >          d = matmul(Hmat,scaled)
864 >
865 >       else
866 >          ! calc the scaled coordinates.
867 >          
868 >          do i = 1, 3
869 >             scaled(i) = d(i) * HmatInv(i,i)
870 >            
871 >             ! wrap the scaled coordinates
872 >            
873 >             scaled(i) = scaled(i) - anint(scaled(i))
874 >            
875 >             ! calc the wrapped real coordinates from the wrapped scaled
876 >             ! coordinates
877 >
878 >             d(i) = scaled(i)*Hmat(i,i)
879 >          enddo
880 >       endif
881        
882      endif
883 +    
884      r_sq = dot_product(d,d)
885 <        
885 >    
886    end subroutine get_interatomic_vector
887 <
887 >  
888    subroutine check_initialization(error)
889      integer, intent(out) :: error
890      
891      error = 0
892      ! Make sure we are properly initialized.
893      if (.not. do_forces_initialized) then
894 +       write(*,*) "Forces not initialized"
895         error = -1
896         return
897      endif
# Line 655 | Line 942 | contains
942      rf = 0.0_dp
943      tau_Temp = 0.0_dp
944      virial_Temp = 0.0_dp
658 –    
945    end subroutine zero_work_arrays
946    
947    function skipThisPair(atom1, atom2) result(skip_it)

Diff Legend

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