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 2512 by gezelter, Thu Dec 15 21:43:16 2005 UTC vs.
Revision 2756 by gezelter, Wed May 17 15:37:15 2006 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.71 2005-12-15 21:43:16 gezelter Exp $, $Date: 2005-12-15 21:43:16 $, $Name: not supported by cvs2svn $, $Revision: 1.71 $
48 > !! @version $Id: doForces.F90,v 1.81 2006-05-17 15:37:14 gezelter Exp $, $Date: 2006-05-17 15:37:14 $, $Name: not supported by cvs2svn $, $Revision: 1.81 $
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 341 | Line 341 | contains
341               endif
342               if (i_is_Shape) then
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
# 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 575 | Line 579 | contains
579          defaultDoShift = .true.
580          
581       endif
582 <
582 >    
583       localError = 0
584       call setLJDefaultCutoff( defaultRcut, defaultDoShift )
585       call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
586 <     call setCutoffEAM( defaultRcut, localError)
587 <     if (localError /= 0) then
588 <       write(errMsg, *) 'An error has occured in setting the EAM cutoff'
589 <       call handleError("setCutoffs", errMsg)
590 <     end if
587 <     call set_switch(GROUP_SWITCH, defaultRsw, defaultRcut)
588 <
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 632 | 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 640 | Line 644 | contains
644  
645    subroutine doReadyCheck(error)
646      integer, intent(out) :: error
643 –
647      integer :: myStatus
648  
649      error = 0
# Line 653 | Line 656 | contains
656         call createGtypeCutoffMap()      
657      endif
658  
656 –
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  
661 –
667      if (.not. haveSIMvariables) then
668         call setSimVariables()
669      endif
670  
666 –  !  if (.not. haveRlist) then
667 –  !     write(default_error, *) 'rList has not been set in doForces!'
668 –  !     error = -1
669 –  !     return
670 –  !  endif
671 –
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 711 | 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 727 | 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 954 | Line 957 | contains
957  
958                     list(nlist) = j
959                  endif
960 <
958 <
959 <                
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, &
969 <                        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 981 | 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 1014 | 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 1175 | 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)
1181 <                   call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1182 <                        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 1238 | Line 1242 | contains
1242      real ( kind = dp ), intent(inout) :: d_grp(3)
1243      real ( kind = dp ), intent(inout) :: rCut
1244      real ( kind = dp ) :: r
1245 +    real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx
1246      integer :: me_i, me_j
1247 +    integer :: k
1248  
1249      integer :: iHash
1250  
1251      r = sqrt(rijsq)
1252 <    vpair = 0.0d0
1253 <    fpair(1:3) = 0.0d0
1252 >    
1253 >    vpair = 0.0_dp
1254 >    fpair(1:3) = 0.0_dp
1255  
1256   #ifdef IS_MPI
1257      me_i = atid_row(i)
# Line 1305 | Line 1312 | contains
1312         call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1313              pot(METALLIC_POT), f, do_pot)
1314      endif
1308 –
1309 –    
1315      
1316    end subroutine do_pair
1317  
# Line 1329 | Line 1334 | contains
1334      integer :: me_i, me_j, iHash
1335  
1336      r = sqrt(rijsq)
1337 <
1337 >    
1338   #ifdef IS_MPI  
1339      me_i = atid_row(i)
1340      me_j = atid_col(j)  
# Line 1361 | Line 1366 | contains
1366      if (FF_uses_SC .and. SIM_uses_SC) then
1367         call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1368      endif
1364 –
1365 –
1369    end subroutine do_preforce
1370  
1371  
# Line 1374 | Line 1377 | contains
1377      real( kind = dp ) :: d(3), scaled(3)
1378      integer i
1379  
1380 <    d(1:3) = q_j(1:3) - q_i(1:3)
1380 >    d(1) = q_j(1) - q_i(1)
1381 >    d(2) = q_j(2) - q_i(2)
1382 >    d(3) = q_j(3) - q_i(3)
1383  
1384      ! Wrap back into periodic box if necessary
1385      if ( SIM_uses_PBC ) then
1386  
1387         if( .not.boxIsOrthorhombic ) then
1388            ! calc the scaled coordinates.
1389 +          ! scaled = matmul(HmatInv, d)
1390  
1391 <          scaled = matmul(HmatInv, d)
1392 <
1391 >          scaled(1) = HmatInv(1,1)*d(1) + HmatInv(1,2)*d(2) + HmatInv(1,3)*d(3)
1392 >          scaled(2) = HmatInv(2,1)*d(1) + HmatInv(2,2)*d(2) + HmatInv(2,3)*d(3)
1393 >          scaled(3) = HmatInv(3,1)*d(1) + HmatInv(3,2)*d(2) + HmatInv(3,3)*d(3)
1394 >          
1395            ! wrap the scaled coordinates
1396  
1397 <          scaled = scaled  - anint(scaled)
1398 <
1397 >          scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1398 >          scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1399 >          scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1400  
1401            ! calc the wrapped real coordinates from the wrapped scaled
1402            ! coordinates
1403 +          ! d = matmul(Hmat,scaled)
1404 +          d(1)= Hmat(1,1)*scaled(1) + Hmat(1,2)*scaled(2) + Hmat(1,3)*scaled(3)
1405 +          d(2)= Hmat(2,1)*scaled(1) + Hmat(2,2)*scaled(2) + Hmat(2,3)*scaled(3)
1406 +          d(3)= Hmat(3,1)*scaled(1) + Hmat(3,2)*scaled(2) + Hmat(3,3)*scaled(3)
1407  
1395 –          d = matmul(Hmat,scaled)
1396 –
1408         else
1409            ! calc the scaled coordinates.
1410  
1411 <          do i = 1, 3
1412 <             scaled(i) = d(i) * HmatInv(i,i)
1411 >          scaled(1) = d(1) * HmatInv(1,1)
1412 >          scaled(2) = d(2) * HmatInv(2,2)
1413 >          scaled(3) = d(3) * HmatInv(3,3)
1414 >          
1415 >          ! wrap the scaled coordinates
1416 >          
1417 >          scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1418 >          scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1419 >          scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1420  
1421 <             ! wrap the scaled coordinates
1421 >          ! calc the wrapped real coordinates from the wrapped scaled
1422 >          ! coordinates
1423  
1424 <             scaled(i) = scaled(i) - anint(scaled(i))
1424 >          d(1) = scaled(1)*Hmat(1,1)
1425 >          d(2) = scaled(2)*Hmat(2,2)
1426 >          d(3) = scaled(3)*Hmat(3,3)
1427  
1407 –             ! calc the wrapped real coordinates from the wrapped scaled
1408 –             ! coordinates
1409 –
1410 –             d(i) = scaled(i)*Hmat(i,i)
1411 –          enddo
1428         endif
1429  
1430      endif
1431  
1432 <    r_sq = dot_product(d,d)
1432 >    r_sq = d(1)*d(1) + d(2)*d(2) + d(3)*d(3)
1433  
1434    end subroutine get_interatomic_vector
1435  

Diff Legend

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