ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/doForces.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2503 by gezelter, Thu Dec 8 22:04:40 2005 UTC vs.
Revision 2787 by gezelter, Mon Jun 5 18:24:45 2006 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.70 2005-12-08 22:04:40 gezelter Exp $, $Date: 2005-12-08 22:04:40 $, $Name: not supported by cvs2svn $, $Revision: 1.70 $
48 > !! @version $Id: doForces.F90,v 1.83 2006-06-05 18:24:45 gezelter Exp $, $Date: 2006-06-05 18:24:45 $, $Name: not supported by cvs2svn $, $Revision: 1.83 $
49  
50  
51   module doForces
# Line 72 | Line 72 | module doForces
72    PRIVATE
73  
74   #define __FORTRAN90
75 – #include "UseTheForce/fSwitchingFunction.h"
75   #include "UseTheForce/fCutoffPolicy.h"
76   #include "UseTheForce/DarkSide/fInteractionMap.h"
77   #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
78  
80 –
79    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
80    INTEGER, PARAMETER:: PAIR_LOOP    = 2
81  
# Line 149 | Line 147 | contains
147    end type gtypeCutoffs
148    type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
149  
150 +
151   contains
152  
153    subroutine createInteractionHash()
# Line 281 | Line 280 | contains
280      logical :: i_is_GB
281      logical :: i_is_EAM
282      logical :: i_is_Shape
283 +    logical :: i_is_SC
284      logical :: GtypeFound
285  
286      integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
# Line 310 | Line 310 | contains
310            call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
311            call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
312            call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
313 <          
313 >          call getElementProperty(atypes, i, "is_SC", i_is_SC)
314  
315            if (haveDefaultCutoffs) then
316               atypeMaxCutoff(i) = defaultRcut
# Line 343 | Line 343 | contains
343                  thisRcut = getShapeCut(i)
344                  if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
345               endif
346 +             if (i_is_SC) then
347 +                thisRcut = getSCCut(i)
348 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
349 +             endif
350            endif
351                      
352            if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
# Line 394 | Line 398 | contains
398         allocate(groupToGtypeCol(jend))
399      end if
400  
401 <    if(.not.associated(groupToGtypeCol)) then
402 <       allocate(groupToGtypeCol(jend))
401 >    if(.not.associated(groupMaxCutoffCol)) then
402 >       allocate(groupMaxCutoffCol(jend))
403      else
404 <       deallocate(groupToGtypeCol)
405 <       allocate(groupToGtypeCol(jend))
404 >       deallocate(groupMaxCutoffCol)
405 >       allocate(groupMaxCutoffCol(jend))
406      end if
407      if(.not.associated(gtypeMaxCutoffCol)) then
408         allocate(gtypeMaxCutoffCol(jend))
# Line 419 | Line 423 | contains
423      !! largest cutoff for any atypes present in this group.  We also
424      !! create gtypes at this point.
425      
426 <    tol = 1.0d-6
426 >    tol = 1.0e-6_dp
427      nGroupTypesRow = 0
428 <
428 >    nGroupTypesCol = 0
429      do i = istart, iend      
430         n_in_i = groupStartRow(i+1) - groupStartRow(i)
431         groupMaxCutoffRow(i) = 0.0_dp
# Line 436 | Line 440 | contains
440               groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
441            endif          
442         enddo
439 –
443         if (nGroupTypesRow.eq.0) then
444            nGroupTypesRow = nGroupTypesRow + 1
445            gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
# Line 576 | Line 579 | contains
579          defaultDoShift = .true.
580          
581       endif
582 <
582 >    
583       localError = 0
584       call setLJDefaultCutoff( defaultRcut, defaultDoShift )
585 <     call setCutoffEAM( defaultRcut, localError)
586 <     if (localError /= 0) then
587 <       write(errMsg, *) 'An error has occured in setting the EAM cutoff'
588 <       call handleError("setCutoffs", errMsg)
589 <     end if
590 <     call set_switch(GROUP_SWITCH, defaultRsw, defaultRcut)
588 <    
585 >     call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
586 >     call setCutoffEAM( defaultRcut )
587 >     call setCutoffSC( defaultRcut )
588 >     call set_switch(defaultRsw, defaultRcut)
589 >     call setHmatDangerousRcutValue(defaultRcut)
590 >        
591       haveDefaultCutoffs = .true.
592 +     haveGtypeCutoffMap = .false.
593 +
594     end subroutine setCutoffs
595  
596     subroutine cWasLame()
# Line 602 | Line 606 | contains
606      
607       cutoffPolicy = cutPolicy
608       haveCutoffPolicy = .true.
609 <     write(*,*) 'have cutoffPolicy in F = ', cutPolicy
606 <
607 <     call createGtypeCutoffMap()
609 >     haveGtypeCutoffMap = .false.
610      
611     end subroutine setCutoffPolicy
612    
# Line 622 | Line 624 | contains
624       real(kind=dp), intent(in) :: thisSkin
625      
626       skinThickness = thisSkin
627 <     haveSkinThickness = .true.
627 >     haveSkinThickness = .true.    
628 >     haveGtypeCutoffMap = .false.
629      
627 –     call createGtypeCutoffMap()
628 –    
630     end subroutine setSkinThickness
631        
632     subroutine setSimVariables()
# Line 634 | Line 635 | contains
635       SIM_requires_postpair_calc = SimRequiresPostpairCalc()
636       SIM_requires_prepair_calc = SimRequiresPrepairCalc()
637       SIM_uses_PBC = SimUsesPBC()
638 +     SIM_uses_SC = SimUsesSC()
639      
640       haveSIMvariables = .true.
641      
# Line 642 | Line 644 | contains
644  
645    subroutine doReadyCheck(error)
646      integer, intent(out) :: error
645 –
647      integer :: myStatus
648  
649      error = 0
# Line 655 | Line 656 | contains
656         call createGtypeCutoffMap()      
657      endif
658  
658 –
659      if (VisitCutoffsAfterComputing) then
660 <       call set_switch(GROUP_SWITCH, largestRcut, largestRcut)      
660 >       call set_switch(largestRcut, largestRcut)      
661 >       call setHmatDangerousRcutValue(largestRcut)
662 >       call setCutoffEAM(largestRcut)
663 >       call setCutoffSC(largestRcut)
664 >       VisitCutoffsAfterComputing = .false.
665      endif
666  
663 –
667      if (.not. haveSIMvariables) then
668         call setSimVariables()
669      endif
670  
668 –  !  if (.not. haveRlist) then
669 –  !     write(default_error, *) 'rList has not been set in doForces!'
670 –  !     error = -1
671 –  !     return
672 –  !  endif
673 –
671      if (.not. haveNeighborList) then
672         write(default_error, *) 'neighbor list has not been initialized in doForces!'
673         error = -1
674         return
675      end if
676 <
676 >    
677      if (.not. haveSaneForceField) then
678         write(default_error, *) 'Force Field is not sane in doForces!'
679         error = -1
680         return
681      end if
682 <
682 >    
683   #ifdef IS_MPI
684      if (.not. isMPISimSet()) then
685         write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
# Line 713 | Line 710 | contains
710      FF_uses_Dipoles = .false.
711      FF_uses_GayBerne = .false.
712      FF_uses_EAM = .false.
713 +    FF_uses_SC = .false.
714  
715      call getMatchingElementList(atypes, "is_Directional", .true., &
716           nMatches, MatchList)
# Line 729 | Line 727 | contains
727      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
728      if (nMatches .gt. 0) FF_uses_EAM = .true.
729  
730 +    call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
731 +    if (nMatches .gt. 0) FF_uses_SC = .true.
732  
733 +
734      haveSaneForceField = .true.
735  
736      if (FF_uses_EAM) then
# Line 956 | Line 957 | contains
957  
958                     list(nlist) = j
959                  endif
960 <
960 <
961 <                
960 >                
961                  if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
962  
963                     rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
964                     if (loop .eq. PAIR_LOOP) then
965 <                      vij = 0.0d0
966 <                      fij(1:3) = 0.0d0
965 >                      vij = 0.0_dp
966 >                      fij(1) = 0.0_dp
967 >                      fij(2) = 0.0_dp
968 >                      fij(3) = 0.0_dp
969                     endif
970                    
971 <                   call get_switch(rgrpsq, sw, dswdr, rgrp, &
971 <                        group_switch, in_switching_region)
971 >                   call get_switch(rgrpsq, sw, dswdr,rgrp, in_switching_region)
972                    
973                     n_in_j = groupStartCol(j+1) - groupStartCol(j)
974                    
# Line 983 | Line 983 | contains
983                           if (skipThisPair(atom1, atom2))  cycle inner
984                          
985                           if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
986 <                            d_atm(1:3) = d_grp(1:3)
986 >                            d_atm(1) = d_grp(1)
987 >                            d_atm(2) = d_grp(2)
988 >                            d_atm(3) = d_grp(3)
989                              ratmsq = rgrpsq
990                           else
991   #ifdef IS_MPI
# Line 1016 | Line 1018 | contains
1018                                   d_grp, rgrp, rCut)
1019   #endif
1020                              vij = vij + vpair
1021 <                            fij(1:3) = fij(1:3) + fpair(1:3)
1021 >                            fij(1) = fij(1) + fpair(1)
1022 >                            fij(2) = fij(2) + fpair(2)
1023 >                            fij(3) = fij(3) + fpair(3)
1024                           endif
1025                        enddo inner
1026                     enddo
# Line 1177 | Line 1181 | contains
1181                  
1182                  ! prevent overcounting of the skips
1183                  if (i.lt.j) then
1184 <                   call get_interatomic_vector(q(:,i), &
1185 <                        q(:,j), d_atm, ratmsq)
1186 <                   rVal = dsqrt(ratmsq)
1183 <                   call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1184 <                        in_switching_region)
1184 >                   call get_interatomic_vector(q(:,i), q(:,j), d_atm, ratmsq)
1185 >                   rVal = sqrt(ratmsq)
1186 >                   call get_switch(ratmsq, sw, dswdr, rVal,in_switching_region)
1187   #ifdef IS_MPI
1188                     call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1189                          vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
# Line 1198 | Line 1200 | contains
1200   #ifdef IS_MPI
1201      
1202      if (do_pot) then
1203 + #ifdef SINGLE_PRECISION
1204 +       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_real,mpi_sum, &
1205 +            mpi_comm_world,mpi_err)            
1206 + #else
1207         call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1208              mpi_comm_world,mpi_err)            
1209 + #endif
1210      endif
1211      
1212      if (do_stress) then
1213 + #ifdef SINGLE_PRECISION
1214 +       call mpi_allreduce(tau_Temp, tau, 9,mpi_real,mpi_sum, &
1215 +            mpi_comm_world,mpi_err)
1216 +       call mpi_allreduce(virial_Temp, virial,1,mpi_real,mpi_sum, &
1217 +            mpi_comm_world,mpi_err)
1218 + #else
1219         call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1220              mpi_comm_world,mpi_err)
1221         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1222              mpi_comm_world,mpi_err)
1223 + #endif
1224      endif
1225      
1226   #else
# Line 1240 | Line 1254 | contains
1254      real ( kind = dp ), intent(inout) :: d_grp(3)
1255      real ( kind = dp ), intent(inout) :: rCut
1256      real ( kind = dp ) :: r
1257 +    real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx
1258      integer :: me_i, me_j
1259 +    integer :: k
1260  
1261      integer :: iHash
1262  
1263      r = sqrt(rijsq)
1264 <    vpair = 0.0d0
1265 <    fpair(1:3) = 0.0d0
1264 >    
1265 >    vpair = 0.0_dp
1266 >    fpair(1:3) = 0.0_dp
1267  
1268   #ifdef IS_MPI
1269      me_i = atid_row(i)
# Line 1284 | Line 1301 | contains
1301      endif
1302      
1303      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1304 <       call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1304 >       call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1305              pot(VDW_POT), A, f, t, do_pot)
1306      endif
1307      
# Line 1307 | Line 1324 | contains
1324         call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1325              pot(METALLIC_POT), f, do_pot)
1326      endif
1310 –
1311 –    
1327      
1328    end subroutine do_pair
1329  
# Line 1331 | Line 1346 | contains
1346      integer :: me_i, me_j, iHash
1347  
1348      r = sqrt(rijsq)
1349 <
1349 >    
1350   #ifdef IS_MPI  
1351      me_i = atid_row(i)
1352      me_j = atid_col(j)  
# Line 1363 | Line 1378 | contains
1378      if (FF_uses_SC .and. SIM_uses_SC) then
1379         call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1380      endif
1366 –
1367 –
1381    end subroutine do_preforce
1382  
1383  
# Line 1376 | Line 1389 | contains
1389      real( kind = dp ) :: d(3), scaled(3)
1390      integer i
1391  
1392 <    d(1:3) = q_j(1:3) - q_i(1:3)
1392 >    d(1) = q_j(1) - q_i(1)
1393 >    d(2) = q_j(2) - q_i(2)
1394 >    d(3) = q_j(3) - q_i(3)
1395  
1396      ! Wrap back into periodic box if necessary
1397      if ( SIM_uses_PBC ) then
1398  
1399         if( .not.boxIsOrthorhombic ) then
1400            ! calc the scaled coordinates.
1401 <
1387 <          scaled = matmul(HmatInv, d)
1401 >          ! scaled = matmul(HmatInv, d)
1402  
1403 +          scaled(1) = HmatInv(1,1)*d(1) + HmatInv(1,2)*d(2) + HmatInv(1,3)*d(3)
1404 +          scaled(2) = HmatInv(2,1)*d(1) + HmatInv(2,2)*d(2) + HmatInv(2,3)*d(3)
1405 +          scaled(3) = HmatInv(3,1)*d(1) + HmatInv(3,2)*d(2) + HmatInv(3,3)*d(3)
1406 +          
1407            ! wrap the scaled coordinates
1408  
1409 <          scaled = scaled  - anint(scaled)
1409 >          scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1410 >          scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1411 >          scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1412  
1393 –
1413            ! calc the wrapped real coordinates from the wrapped scaled
1414            ! coordinates
1415 +          ! d = matmul(Hmat,scaled)
1416 +          d(1)= Hmat(1,1)*scaled(1) + Hmat(1,2)*scaled(2) + Hmat(1,3)*scaled(3)
1417 +          d(2)= Hmat(2,1)*scaled(1) + Hmat(2,2)*scaled(2) + Hmat(2,3)*scaled(3)
1418 +          d(3)= Hmat(3,1)*scaled(1) + Hmat(3,2)*scaled(2) + Hmat(3,3)*scaled(3)
1419  
1397 –          d = matmul(Hmat,scaled)
1398 –
1420         else
1421            ! calc the scaled coordinates.
1422  
1423 <          do i = 1, 3
1424 <             scaled(i) = d(i) * HmatInv(i,i)
1423 >          scaled(1) = d(1) * HmatInv(1,1)
1424 >          scaled(2) = d(2) * HmatInv(2,2)
1425 >          scaled(3) = d(3) * HmatInv(3,3)
1426 >          
1427 >          ! wrap the scaled coordinates
1428 >          
1429 >          scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1430 >          scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1431 >          scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1432  
1433 <             ! wrap the scaled coordinates
1433 >          ! calc the wrapped real coordinates from the wrapped scaled
1434 >          ! coordinates
1435  
1436 <             scaled(i) = scaled(i) - anint(scaled(i))
1436 >          d(1) = scaled(1)*Hmat(1,1)
1437 >          d(2) = scaled(2)*Hmat(2,2)
1438 >          d(3) = scaled(3)*Hmat(3,3)
1439  
1409 –             ! calc the wrapped real coordinates from the wrapped scaled
1410 –             ! coordinates
1411 –
1412 –             d(i) = scaled(i)*Hmat(i,i)
1413 –          enddo
1440         endif
1441  
1442      endif
1443  
1444 <    r_sq = dot_product(d,d)
1444 >    r_sq = d(1)*d(1) + d(2)*d(2) + d(3)*d(3)
1445  
1446    end subroutine get_interatomic_vector
1447  

Diff Legend

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