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 2530 by chuckv, Fri Dec 30 00:18:28 2005 UTC vs.
Revision 2722 by gezelter, Thu Apr 20 18:24:24 2006 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.72 2005-12-30 00:18:28 chuckv Exp $, $Date: 2005-12-30 00:18:28 $, $Name: not supported by cvs2svn $, $Revision: 1.72 $
48 > !! @version $Id: doForces.F90,v 1.79 2006-04-20 18:24:24 gezelter Exp $, $Date: 2006-04-20 18:24:24 $, $Name: not supported by cvs2svn $, $Revision: 1.79 $
49  
50  
51   module doForces
# Line 64 | Line 64 | module doForces
64    use eam
65    use suttonchen
66    use status
67 +  use interpolation
68   #ifdef IS_MPI
69    use mpiSimulation
70   #endif
# Line 72 | Line 73 | module doForces
73    PRIVATE
74  
75   #define __FORTRAN90
75 – #include "UseTheForce/fSwitchingFunction.h"
76   #include "UseTheForce/fCutoffPolicy.h"
77   #include "UseTheForce/DarkSide/fInteractionMap.h"
78   #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
79  
80 –
80    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
81    INTEGER, PARAMETER:: PAIR_LOOP    = 2
82 +  INTEGER, PARAMETER:: np = 500
83  
84    logical, save :: haveNeighborList = .false.
85    logical, save :: haveSIMvariables = .false.
# Line 91 | Line 91 | module doForces
91    logical, save :: haveElectrostaticSummationMethod = .false.
92    logical, save :: haveCutoffPolicy = .false.
93    logical, save :: VisitCutoffsAfterComputing = .false.
94 +  logical, save :: haveSplineSqrt = .false.
95  
96    logical, save :: FF_uses_DirectionalAtoms
97    logical, save :: FF_uses_Dipoles
# Line 149 | Line 150 | contains
150    end type gtypeCutoffs
151    type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
152  
153 +  ! variables for the spline of the sqrt
154 +  type(cubicSpline), save :: splineSqrt
155 +  logical, save :: useSpline = .true.
156 +  
157 +
158   contains
159  
160    subroutine createInteractionHash()
# Line 399 | Line 405 | contains
405         allocate(groupToGtypeCol(jend))
406      end if
407  
408 <    if(.not.associated(groupToGtypeCol)) then
409 <       allocate(groupToGtypeCol(jend))
408 >    if(.not.associated(groupMaxCutoffCol)) then
409 >       allocate(groupMaxCutoffCol(jend))
410      else
411 <       deallocate(groupToGtypeCol)
412 <       allocate(groupToGtypeCol(jend))
411 >       deallocate(groupMaxCutoffCol)
412 >       allocate(groupMaxCutoffCol(jend))
413      end if
414      if(.not.associated(gtypeMaxCutoffCol)) then
415         allocate(gtypeMaxCutoffCol(jend))
# Line 426 | Line 432 | contains
432      
433      tol = 1.0d-6
434      nGroupTypesRow = 0
435 <
435 >    nGroupTypesCol = 0
436      do i = istart, iend      
437         n_in_i = groupStartRow(i+1) - groupStartRow(i)
438         groupMaxCutoffRow(i) = 0.0_dp
# Line 580 | Line 586 | contains
586          defaultDoShift = .true.
587          
588       endif
589 <
589 >    
590       localError = 0
591       call setLJDefaultCutoff( defaultRcut, defaultDoShift )
592       call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
593 <     call setCutoffEAM( defaultRcut, localError)
594 <     if (localError /= 0) then
595 <       write(errMsg, *) 'An error has occured in setting the EAM cutoff'
596 <       call handleError("setCutoffs", errMsg)
597 <     end if
598 <     call set_switch(GROUP_SWITCH, defaultRsw, defaultRcut)
593 <
593 >     call setCutoffEAM( defaultRcut )
594 >     call setCutoffSC( defaultRcut )
595 >     call set_switch(defaultRsw, defaultRcut)
596 >     call setHmatDangerousRcutValue(defaultRcut)
597 >     call setupSplineSqrt(defaultRcut)
598 >        
599       haveDefaultCutoffs = .true.
600       haveGtypeCutoffMap = .false.
601 +
602     end subroutine setCutoffs
603  
604     subroutine cWasLame()
# Line 637 | Line 643 | contains
643       SIM_requires_postpair_calc = SimRequiresPostpairCalc()
644       SIM_requires_prepair_calc = SimRequiresPrepairCalc()
645       SIM_uses_PBC = SimUsesPBC()
646 +     SIM_uses_SC = SimUsesSC()
647      
648       haveSIMvariables = .true.
649      
# Line 645 | Line 652 | contains
652  
653    subroutine doReadyCheck(error)
654      integer, intent(out) :: error
648 –
655      integer :: myStatus
656  
657      error = 0
# Line 658 | Line 664 | contains
664         call createGtypeCutoffMap()      
665      endif
666  
661 –
667      if (VisitCutoffsAfterComputing) then
668 <       call set_switch(GROUP_SWITCH, largestRcut, largestRcut)      
668 >       call set_switch(largestRcut, largestRcut)      
669 >       call setHmatDangerousRcutValue(largestRcut)
670 >       call setLJsplineRmax(largestRcut)
671 >       call setCutoffEAM(largestRcut)
672 >       call setCutoffSC(largestRcut)
673 >       VisitCutoffsAfterComputing = .false.
674      endif
675  
676 +    if (.not. haveSplineSqrt) then
677 +       call setupSplineSqrt(largestRcut)
678 +    endif
679  
680      if (.not. haveSIMvariables) then
681         call setSimVariables()
682      endif
683  
671 –  !  if (.not. haveRlist) then
672 –  !     write(default_error, *) 'rList has not been set in doForces!'
673 –  !     error = -1
674 –  !     return
675 –  !  endif
676 –
684      if (.not. haveNeighborList) then
685         write(default_error, *) 'neighbor list has not been initialized in doForces!'
686         error = -1
687         return
688      end if
689 <
689 >    
690      if (.not. haveSaneForceField) then
691         write(default_error, *) 'Force Field is not sane in doForces!'
692         error = -1
693         return
694      end if
695 <
695 >    
696   #ifdef IS_MPI
697      if (.not. isMPISimSet()) then
698         write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
# Line 716 | Line 723 | contains
723      FF_uses_Dipoles = .false.
724      FF_uses_GayBerne = .false.
725      FF_uses_EAM = .false.
726 +    FF_uses_SC = .false.
727  
728      call getMatchingElementList(atypes, "is_Directional", .true., &
729           nMatches, MatchList)
# Line 732 | Line 740 | contains
740      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
741      if (nMatches .gt. 0) FF_uses_EAM = .true.
742  
743 +    call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
744 +    if (nMatches .gt. 0) FF_uses_SC = .true.
745  
746 +
747      haveSaneForceField = .true.
748  
749      if (FF_uses_EAM) then
# Line 959 | Line 970 | contains
970  
971                     list(nlist) = j
972                  endif
973 <
963 <
964 <                
973 >                
974                  if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
975  
976                     rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
977                     if (loop .eq. PAIR_LOOP) then
978                        vij = 0.0d0
979 <                      fij(1:3) = 0.0d0
979 >                      fij(1) = 0.0_dp
980 >                      fij(2) = 0.0_dp
981 >                      fij(3) = 0.0_dp
982                     endif
983                    
984 <                   call get_switch(rgrpsq, sw, dswdr, rgrp, &
974 <                        group_switch, in_switching_region)
984 >                   call get_switch(rgrpsq, sw, dswdr,rgrp, in_switching_region)
985                    
986                     n_in_j = groupStartCol(j+1) - groupStartCol(j)
987                    
# Line 986 | Line 996 | contains
996                           if (skipThisPair(atom1, atom2))  cycle inner
997                          
998                           if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
999 <                            d_atm(1:3) = d_grp(1:3)
999 >                            d_atm(1) = d_grp(1)
1000 >                            d_atm(2) = d_grp(2)
1001 >                            d_atm(3) = d_grp(3)
1002                              ratmsq = rgrpsq
1003                           else
1004   #ifdef IS_MPI
# Line 1019 | Line 1031 | contains
1031                                   d_grp, rgrp, rCut)
1032   #endif
1033                              vij = vij + vpair
1034 <                            fij(1:3) = fij(1:3) + fpair(1:3)
1034 >                            fij(1) = fij(1) + fpair(1)
1035 >                            fij(2) = fij(2) + fpair(2)
1036 >                            fij(3) = fij(3) + fpair(3)
1037                           endif
1038                        enddo inner
1039                     enddo
# Line 1180 | Line 1194 | contains
1194                  
1195                  ! prevent overcounting of the skips
1196                  if (i.lt.j) then
1197 <                   call get_interatomic_vector(q(:,i), &
1184 <                        q(:,j), d_atm, ratmsq)
1197 >                   call get_interatomic_vector(q(:,i), q(:,j), d_atm, ratmsq)
1198                     rVal = dsqrt(ratmsq)
1199 <                   call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1187 <                        in_switching_region)
1199 >                   call get_switch(ratmsq, sw, dswdr, rVal,in_switching_region)
1200   #ifdef IS_MPI
1201                     call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1202                          vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
# Line 1243 | Line 1255 | contains
1255      real ( kind = dp ), intent(inout) :: d_grp(3)
1256      real ( kind = dp ), intent(inout) :: rCut
1257      real ( kind = dp ) :: r
1258 +    real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx
1259      integer :: me_i, me_j
1260 +    integer :: k
1261  
1262      integer :: iHash
1263  
1264 <    r = sqrt(rijsq)
1264 >    if (useSpline) then
1265 >       call lookupUniformSpline(splineSqrt, rijsq, r)
1266 >    else
1267 >       r = sqrt(rijsq)
1268 >    endif
1269 >
1270      vpair = 0.0d0
1271      fpair(1:3) = 0.0d0
1272  
# Line 1310 | Line 1329 | contains
1329         call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1330              pot(METALLIC_POT), f, do_pot)
1331      endif
1313 –
1314 –    
1332      
1333    end subroutine do_pair
1334  
# Line 1333 | Line 1350 | contains
1350  
1351      integer :: me_i, me_j, iHash
1352  
1353 <    r = sqrt(rijsq)
1353 >    if (useSpline) then
1354 >       call lookupUniformSpline(splineSqrt, rijsq, r)
1355 >    else
1356 >       r = sqrt(rijsq)
1357 >    endif
1358  
1359   #ifdef IS_MPI  
1360      me_i = atid_row(i)
# Line 1366 | Line 1387 | contains
1387      if (FF_uses_SC .and. SIM_uses_SC) then
1388         call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1389      endif
1369 –
1370 –
1390    end subroutine do_preforce
1391  
1392  
# Line 1379 | Line 1398 | contains
1398      real( kind = dp ) :: d(3), scaled(3)
1399      integer i
1400  
1401 <    d(1:3) = q_j(1:3) - q_i(1:3)
1401 >    d(1) = q_j(1) - q_i(1)
1402 >    d(2) = q_j(2) - q_i(2)
1403 >    d(3) = q_j(3) - q_i(3)
1404  
1405      ! Wrap back into periodic box if necessary
1406      if ( SIM_uses_PBC ) then
1407  
1408         if( .not.boxIsOrthorhombic ) then
1409            ! calc the scaled coordinates.
1410 +          ! scaled = matmul(HmatInv, d)
1411  
1412 <          scaled = matmul(HmatInv, d)
1413 <
1412 >          scaled(1) = HmatInv(1,1)*d(1) + HmatInv(1,2)*d(2) + HmatInv(1,3)*d(3)
1413 >          scaled(2) = HmatInv(2,1)*d(1) + HmatInv(2,2)*d(2) + HmatInv(2,3)*d(3)
1414 >          scaled(3) = HmatInv(3,1)*d(1) + HmatInv(3,2)*d(2) + HmatInv(3,3)*d(3)
1415 >          
1416            ! wrap the scaled coordinates
1417  
1418 <          scaled = scaled  - anint(scaled)
1418 >          scaled(1) = scaled(1) - dnint(scaled(1))
1419 >          scaled(2) = scaled(2) - dnint(scaled(2))
1420 >          scaled(3) = scaled(3) - dnint(scaled(3))
1421  
1396 –
1422            ! calc the wrapped real coordinates from the wrapped scaled
1423            ! coordinates
1424 +          ! d = matmul(Hmat,scaled)
1425 +          d(1)= Hmat(1,1)*scaled(1) + Hmat(1,2)*scaled(2) + Hmat(1,3)*scaled(3)
1426 +          d(2)= Hmat(2,1)*scaled(1) + Hmat(2,2)*scaled(2) + Hmat(2,3)*scaled(3)
1427 +          d(3)= Hmat(3,1)*scaled(1) + Hmat(3,2)*scaled(2) + Hmat(3,3)*scaled(3)
1428  
1400 –          d = matmul(Hmat,scaled)
1401 –
1429         else
1430            ! calc the scaled coordinates.
1431  
1432 <          do i = 1, 3
1433 <             scaled(i) = d(i) * HmatInv(i,i)
1434 <
1435 <             ! wrap the scaled coordinates
1436 <
1437 <             scaled(i) = scaled(i) - anint(scaled(i))
1432 >          scaled(1) = d(1) * HmatInv(1,1)
1433 >          scaled(2) = d(2) * HmatInv(2,2)
1434 >          scaled(3) = d(3) * HmatInv(3,3)
1435 >          
1436 >          ! wrap the scaled coordinates
1437 >          
1438 >          scaled(1) = scaled(1) - dnint(scaled(1))
1439 >          scaled(2) = scaled(2) - dnint(scaled(2))
1440 >          scaled(3) = scaled(3) - dnint(scaled(3))
1441  
1442 <             ! calc the wrapped real coordinates from the wrapped scaled
1443 <             ! coordinates
1442 >          ! calc the wrapped real coordinates from the wrapped scaled
1443 >          ! coordinates
1444  
1445 <             d(i) = scaled(i)*Hmat(i,i)
1446 <          enddo
1445 >          d(1) = scaled(1)*Hmat(1,1)
1446 >          d(2) = scaled(2)*Hmat(2,2)
1447 >          d(3) = scaled(3)*Hmat(3,3)
1448 >
1449         endif
1450  
1451      endif
1452  
1453 <    r_sq = dot_product(d,d)
1453 >    r_sq = d(1)*d(1) + d(2)*d(2) + d(3)*d(3)
1454  
1455    end subroutine get_interatomic_vector
1456  
# Line 1585 | Line 1617 | end module doForces
1617  
1618    end subroutine add_stress_tensor
1619  
1620 +  subroutine setupSplineSqrt(rmax)
1621 +    real(kind=dp), intent(in) :: rmax
1622 +    real(kind=dp), dimension(np) :: xvals, yvals
1623 +    real(kind=dp) :: r2_1, r2_n, dx, r2
1624 +    integer :: i
1625 +
1626 +    r2_1 = 0.5d0
1627 +    r2_n = rmax*rmax
1628 +
1629 +    dx = (r2_n-r2_1) / dble(np-1)
1630 +    
1631 +    do i = 1, np
1632 +       r2 = r2_1 + dble(i-1)*dx
1633 +       xvals(i) = r2
1634 +       yvals(i) = dsqrt(r2)
1635 +    enddo
1636 +
1637 +    call newSpline(splineSqrt, xvals, yvals, .true.)
1638 +    
1639 +    haveSplineSqrt = .true.
1640 +    return
1641 +  end subroutine setupSplineSqrt
1642 +
1643 +  subroutine deleteSplineSqrt()
1644 +    call deleteSpline(splineSqrt)
1645 +    haveSplineSqrt = .false.
1646 +    return
1647 +  end subroutine deleteSplineSqrt
1648 +
1649   end module doForces

Diff Legend

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