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 486 by mmeineke, Thu Apr 10 16:22:00 2003 UTC vs.
Revision 650 by chuckv, Thu Jul 24 19:57:35 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.14 2003-04-10 16:22:00 mmeineke Exp $, $Date: 2003-04-10 16:22:00 $, $Name: not supported by cvs2svn $, $Revision: 1.14 $
7 > !! @version $Id: do_Forces.F90,v 1.25 2003-07-24 19:57:35 chuckv Exp $, $Date: 2003-07-24 19:57:35 $, $Name: not supported by cvs2svn $, $Revision: 1.25 $
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   #ifdef IS_MPI
23    use mpiSimulation
24   #endif
# Line 27 | Line 29 | module do_Forces
29   #define __FORTRAN90
30   #include "fForceField.h"
31  
32 <  logical, save :: do_forces_initialized = .false.
32 >  logical, save :: do_forces_initialized = .false., haveRlist = .false.
33 >  logical, save :: havePolicies = .false.
34    logical, save :: FF_uses_LJ
35    logical, save :: FF_uses_sticky
36    logical, save :: FF_uses_dipoles
# Line 35 | Line 38 | module do_Forces
38    logical, save :: FF_uses_GB
39    logical, save :: FF_uses_EAM
40  
41 +  real(kind=dp), save :: rlist, rlistsq
42 +
43    public :: init_FF
44    public :: do_force_loop
45 +  public :: setRlistDF
46  
47   contains
48  
49 +  subroutine setRlistDF( this_rlist )
50 +    
51 +    real(kind=dp) :: this_rlist
52 +
53 +    rlist = this_rlist
54 +    rlistsq = rlist * rlist
55 +    
56 +    haveRlist = .true.
57 +    if( havePolicies ) do_forces_initialized = .true.
58 +
59 +  end subroutine setRlistDF    
60 +
61    subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
62  
63      integer, intent(in) :: LJMIXPOLICY
# Line 87 | Line 105 | contains
105      !! check to make sure the FF_uses_RF setting makes sense
106      
107      if (FF_uses_dipoles) then
90 –       rrf = getRrf()
91 –       rt = getRt()      
92 –       call initialize_dipole(rrf, rt)
108         if (FF_uses_RF) then
109            dielect = getDielect()
110 <          call initialize_rf(rrf, rt, dielect)
110 >          call initialize_rf(dielect)
111         endif
112      else
113         if (FF_uses_RF) then          
# Line 100 | Line 115 | contains
115            thisStat = -1
116            return
117         endif
118 <    endif
118 >    endif
119  
120      if (FF_uses_LJ) then
121        
107 –       call getRcut(rcut)
108 –
122         select case (LJMIXPOLICY)
123         case (LB_MIXING_RULE)
124 <          call init_lj_FF(LB_MIXING_RULE, rcut, my_status)            
124 >          call init_lj_FF(LB_MIXING_RULE, my_status)            
125         case (EXPLICIT_MIXING_RULE)
126 <          call init_lj_FF(EXPLICIT_MIXING_RULE, rcut, my_status)
126 >          call init_lj_FF(EXPLICIT_MIXING_RULE, my_status)
127         case default
128            write(default_error,*) 'unknown LJ Mixing Policy!'
129            thisStat = -1
# Line 150 | Line 163 | contains
163         endif
164      endif
165  
166 <    do_forces_initialized = .true.    
167 <
166 >    havePolicies = .true.
167 >    if( haveRlist ) do_forces_initialized = .true.
168 >    
169    end subroutine init_FF
170    
171  
# Line 185 | Line 199 | contains
199      logical :: update_nlist  
200      integer :: i, j, jbeg, jend, jnab
201      integer :: nlist
202 <    real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut
202 >    real( kind = DP ) ::  rijsq
203      real(kind=dp),dimension(3) :: d
204      real(kind=dp) :: rfpot, mu_i, virial
205      integer :: me_i
# Line 193 | Line 207 | contains
207      integer :: neighborListSize
208      integer :: listerror, error
209      integer :: localError
210 +
211 +    real(kind=dp) :: listSkin = 1.0
212 +    
213  
214      !! initialize local variables  
215  
# Line 206 | Line 223 | contains
223      natoms = nlocal
224   #endif
225    
209 –    call getRcut(rcut,rc2=rcutsq)
210 –    call getRlist(rlist,rlistsq)
211 –    
226      call check_initialization(localError)
227      if ( localError .ne. 0 ) then
228         error = -1
# Line 238 | Line 252 | contains
252      
253      if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
254         !! See if we need to update neighbor lists
255 <       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)  
255 >       call checkNeighborList(nlocal, q, listSkin, update_nlist)  
256         !! if_mpi_gather_stuff_for_prepair
257         !! do_prepair_loop_if_needed
258         !! if_mpi_scatter_stuff_from_prepair
259         !! if_mpi_gather_stuff_from_prepair_to_main_loop
260 <    else
261 <       !! See if we need to update neighbor lists
262 <       call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)  
260 >
261 > !--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
262 > #ifdef IS_MPI
263 >    
264 >    if (update_nlist) then
265 >      
266 >       !! save current configuration, construct neighbor list,
267 >       !! and calculate forces
268 >       call saveNeighborList(nlocal, q)
269 >      
270 >       neighborListSize = size(list)
271 >       nlist = 0      
272 >      
273 >       do i = 1, nrow
274 >          point(i) = nlist + 1
275 >          
276 >          prepair_inner: do j = 1, ncol
277 >            
278 >             if (skipThisPair(i,j)) cycle prepair_inner
279 >            
280 >             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
281 >            
282 >             if (rijsq < rlistsq) then            
283 >                
284 >                nlist = nlist + 1
285 >                
286 >                if (nlist > neighborListSize) then
287 >                   call expandNeighborList(nlocal, listerror)
288 >                   if (listerror /= 0) then
289 >                      error = -1
290 >                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
291 >                      return
292 >                   end if
293 >                   neighborListSize = size(list)
294 >                endif
295 >                
296 >                list(nlist) = j
297 >                call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)                      
298 >             endif
299 >          enddo prepair_inner
300 >       enddo
301 >
302 >       point(nrow + 1) = nlist + 1
303 >      
304 >    else  !! (of update_check)
305 >
306 >       ! use the list to find the neighbors
307 >       do i = 1, nrow
308 >          JBEG = POINT(i)
309 >          JEND = POINT(i+1) - 1
310 >          ! check thiat molecule i has neighbors
311 >          if (jbeg .le. jend) then
312 >            
313 >             do jnab = jbeg, jend
314 >                j = list(jnab)
315 >
316 >                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
317 >                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
318 >                     u_l, A, f, t, pot_local)
319 >
320 >             enddo
321 >          endif
322 >       enddo
323      endif
324      
325 + #else
326 +    
327 +    if (update_nlist) then
328 +      
329 +       ! save current configuration, contruct neighbor list,
330 +       ! and calculate forces
331 +       call saveNeighborList(natoms, q)
332 +      
333 +       neighborListSize = size(list)
334 +  
335 +       nlist = 0
336 +      
337 +       do i = 1, natoms-1
338 +          point(i) = nlist + 1
339 +          
340 +          prepair_inner: do j = i+1, natoms
341 +            
342 +             if (skipThisPair(i,j))  cycle prepair_inner
343 +                          
344 +             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
345 +          
346 +
347 +             if (rijsq < rlistsq) then
348 +                
349 +                nlist = nlist + 1
350 +              
351 +                if (nlist > neighborListSize) then
352 +                   call expandNeighborList(natoms, listerror)
353 +                   if (listerror /= 0) then
354 +                      error = -1
355 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
356 +                      return
357 +                   end if
358 +                   neighborListSize = size(list)
359 +                endif
360 +                
361 +                list(nlist) = j
362 +                
363 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
364 +                        u_l, A, f, t, pot)
365 +                
366 +             endif
367 +          enddo prepair_inner
368 +       enddo
369 +      
370 +       point(natoms) = nlist + 1
371 +      
372 +    else !! (update)
373 +      
374 +       ! use the list to find the neighbors
375 +       do i = 1, natoms-1
376 +          JBEG = POINT(i)
377 +          JEND = POINT(i+1) - 1
378 +          ! check thiat molecule i has neighbors
379 +          if (jbeg .le. jend) then
380 +            
381 +             do jnab = jbeg, jend
382 +                j = list(jnab)
383 +
384 +                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
385 +                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
386 +                     u_l, A, f, t, pot)
387 +
388 +             enddo
389 +          endif
390 +       enddo
391 +    endif    
392 + #endif
393 +    !! Do rest of preforce calculations
394 +   call do_preforce(nlocal,pot)
395 +    else
396 +       !! See if we need to update neighbor lists for non pre-pair
397 +       call checkNeighborList(nlocal, q, listSkin, update_nlist)  
398 +    endif
399 +
400 +
401 +
402 +
403 +
404 + !---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
405 +
406 +
407 +
408 +
409 +  
410   #ifdef IS_MPI
411      
412      if (update_nlist) then
# Line 268 | Line 427 | contains
427              
428               call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
429              
430 <             if (rijsq <  rlistsq) then            
430 >             if (rijsq < rlistsq) then            
431                  
432                  nlist = nlist + 1
433                  
# Line 284 | Line 443 | contains
443                  
444                  list(nlist) = j
445                                  
446 <                if (rijsq <  rcutsq) then
447 <                   call do_pair(i, j, rijsq, d, do_pot, do_stress, &
448 <                        u_l, A, f, t, pot_local)
290 <                endif
446 >                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
447 >                     u_l, A, f, t, pot_local)
448 >                
449               endif
450            enddo inner
451         enddo
# Line 337 | Line 495 | contains
495               call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
496            
497  
498 <             if (rijsq <  rlistsq) then
498 >             if (rijsq < rlistsq) then
499                  
500                  nlist = nlist + 1
501                
# Line 353 | Line 511 | contains
511                  
512                  list(nlist) = j
513                  
514 <                if (rijsq <  rcutsq) then
357 <                   call do_pair(i, j, rijsq, d, do_pot, do_stress, &
514 >                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
515                          u_l, A, f, t, pot)
516 <                endif
516 >                
517               endif
518            enddo inner
519         enddo
# Line 487 | Line 644 | contains
644      endif
645  
646      if (do_stress) then
647 <       call mpi_allreduce(tau_Temp, tau,9,mpi_double_precision,mpi_sum, &
647 >      call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
648              mpi_comm_world,mpi_err)
649         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
650              mpi_comm_world,mpi_err)
# Line 520 | Line 677 | contains
677      logical :: is_LJ_i, is_LJ_j
678      logical :: is_DP_i, is_DP_j
679      logical :: is_GB_i, is_GB_j
680 +    logical :: is_EAM_i,is_EAM_j
681      logical :: is_Sticky_i, is_Sticky_j
682      integer :: me_i, me_j
683  
684      r = sqrt(rijsq)
685  
686   #ifdef IS_MPI
687 +    if (tagRow(i) .eq. tagColumn(j)) then
688 +       write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
689 +    endif
690  
691      me_i = atid_row(i)
692      me_j = atid_col(j)
# Line 584 | Line 745 | contains
745         endif
746      endif
747      
748 +
749 +  
750 +   if (FF_uses_EAM .and. SimUsesEAM()) then
751 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
752 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
753 +      
754 +      if ( is_EAM_i .and. is_EAM_j ) &
755 +           call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
756 +   endif
757 +
758 +
759 +
760 +
761    end subroutine do_pair
762 +
763 +
764 +
765 +  subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
766 +   real( kind = dp ) :: pot
767 +   real( kind = dp ), dimension(3,getNlocal()) :: u_l
768 +   real (kind=dp), dimension(9,getNlocal()) :: A
769 +   real (kind=dp), dimension(3,getNlocal()) :: f
770 +   real (kind=dp), dimension(3,getNlocal()) :: t
771 +  
772 +   logical, intent(inout) :: do_pot, do_stress
773 +   integer, intent(in) :: i, j
774 +   real ( kind = dp ), intent(inout)    :: rijsq
775 +   real ( kind = dp )                :: r
776 +   real ( kind = dp ), intent(inout) :: d(3)
777 +  
778 +   logical :: is_EAM_i, is_EAM_j
779 +  
780 +   integer :: me_i, me_j
781 +  
782 +   r = sqrt(rijsq)
783 +  
784 + #ifdef IS_MPI
785 +   if (tagRow(i) .eq. tagColumn(j)) then
786 +      write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
787 +   endif
788 +  
789 +   me_i = atid_row(i)
790 +   me_j = atid_col(j)
791 +  
792 + #else
793 +  
794 +   me_i = atid(i)
795 +   me_j = atid(j)
796 +  
797 + #endif
798 +  
799 +   if (FF_uses_EAM .and. SimUsesEAM()) then
800 +      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
801 +      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
802 +      
803 +      if ( is_EAM_i .and. is_EAM_j ) &
804 +           call calc_EAM_prepair_rho(i, j, d, r, rijsq )
805 +   endif
806 +  end subroutine do_prepair
807  
808  
809 +
810 +
811 +  subroutine do_preforce(nlocal,pot)
812 +    integer :: nlocal
813 +    real( kind = dp ) :: pot
814 +
815 +   if (FF_uses_EAM .and. SimUsesEAM()) then
816 +      call calc_EAM_preforce_Frho(nlocal,pot)
817 +   endif
818 +
819 +
820 +  end subroutine do_preforce
821 +  
822 +  
823    subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
824      
825      real (kind = dp), dimension(3) :: q_i
826      real (kind = dp), dimension(3) :: q_j
827      real ( kind = dp ), intent(out) :: r_sq
828 <    real( kind = dp ) :: d(3)
829 <    real( kind = dp ) :: d_old(3)
828 >    real( kind = dp ) :: d(3), scaled(3)
829 >    integer i
830 >
831      d(1:3) = q_j(1:3) - q_i(1:3)
832 <    d_old = d
832 >
833      ! Wrap back into periodic box if necessary
834      if ( SimUsesPBC() ) then
835        
836 <       d(1:3) = d(1:3) - box(1:3) * sign(1.0_dp,d(1:3)) * &
837 <            int(abs(d(1:3)/box(1:3)) + 0.5_dp)
836 >       if( .not.boxIsOrthorhombic ) then
837 >          ! calc the scaled coordinates.
838 >          
839 >          scaled = matmul(HmatInv, d)
840 >          
841 >          ! wrap the scaled coordinates
842 >
843 >          scaled = scaled  - anint(scaled)
844 >          
845 >
846 >          ! calc the wrapped real coordinates from the wrapped scaled
847 >          ! coordinates
848 >
849 >          d = matmul(Hmat,scaled)
850 >
851 >       else
852 >          ! calc the scaled coordinates.
853 >          
854 >          do i = 1, 3
855 >             scaled(i) = d(i) * HmatInv(i,i)
856 >            
857 >             ! wrap the scaled coordinates
858 >            
859 >             scaled(i) = scaled(i) - anint(scaled(i))
860 >            
861 >             ! calc the wrapped real coordinates from the wrapped scaled
862 >             ! coordinates
863 >
864 >             d(i) = scaled(i)*Hmat(i,i)
865 >          enddo
866 >       endif
867        
868      endif
869 +    
870      r_sq = dot_product(d,d)
871 <        
871 >    
872    end subroutine get_interatomic_vector
873 <
873 >  
874    subroutine check_initialization(error)
875      integer, intent(out) :: error
876      

Diff Legend

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