| 45 |
|
|
| 46 |
|
!! @author Charles F. Vardeman II |
| 47 |
|
!! @author Matthew Meineke |
| 48 |
< |
!! @version $Id: doForces.F90,v 1.73 2005-12-30 21:25:56 tim Exp $, $Date: 2005-12-30 21:25:56 $, $Name: not supported by cvs2svn $, $Revision: 1.73 $ |
| 48 |
> |
!! @version $Id: doForces.F90,v 1.78 2006-04-17 21:49:12 gezelter Exp $, $Date: 2006-04-17 21:49:12 $, $Name: not supported by cvs2svn $, $Revision: 1.78 $ |
| 49 |
|
|
| 50 |
|
|
| 51 |
|
module doForces |
| 64 |
|
use eam |
| 65 |
|
use suttonchen |
| 66 |
|
use status |
| 67 |
+ |
use interpolation |
| 68 |
|
#ifdef IS_MPI |
| 69 |
|
use mpiSimulation |
| 70 |
|
#endif |
| 585 |
|
localError = 0 |
| 586 |
|
call setLJDefaultCutoff( defaultRcut, defaultDoShift ) |
| 587 |
|
call setElectrostaticCutoffRadius( defaultRcut, defaultRsw ) |
| 588 |
< |
call setCutoffEAM( defaultRcut, localError) |
| 589 |
< |
if (localError /= 0) then |
| 589 |
< |
write(errMsg, *) 'An error has occured in setting the EAM cutoff' |
| 590 |
< |
call handleError("setCutoffs", errMsg) |
| 591 |
< |
end if |
| 588 |
> |
call setCutoffEAM( defaultRcut ) |
| 589 |
> |
call setCutoffSC( defaultRcut ) |
| 590 |
|
call set_switch(GROUP_SWITCH, defaultRsw, defaultRcut) |
| 591 |
+ |
call setHmatDangerousRcutValue(defaultRcut) |
| 592 |
|
|
| 593 |
|
haveDefaultCutoffs = .true. |
| 594 |
|
haveGtypeCutoffMap = .false. |
| 636 |
|
SIM_requires_postpair_calc = SimRequiresPostpairCalc() |
| 637 |
|
SIM_requires_prepair_calc = SimRequiresPrepairCalc() |
| 638 |
|
SIM_uses_PBC = SimUsesPBC() |
| 639 |
+ |
SIM_uses_SC = SimUsesSC() |
| 640 |
|
|
| 641 |
|
haveSIMvariables = .true. |
| 642 |
|
|
| 661 |
|
|
| 662 |
|
if (VisitCutoffsAfterComputing) then |
| 663 |
|
call set_switch(GROUP_SWITCH, largestRcut, largestRcut) |
| 664 |
+ |
call setHmatDangerousRcutValue(largestRcut) |
| 665 |
+ |
call setLJsplineRmax(largestRcut) |
| 666 |
+ |
call setCutoffEAM(largestRcut) |
| 667 |
+ |
call setCutoffSC(largestRcut) |
| 668 |
+ |
VisitCutoffsAfterComputing = .false. |
| 669 |
|
endif |
| 670 |
|
|
| 671 |
|
|
| 721 |
|
FF_uses_Dipoles = .false. |
| 722 |
|
FF_uses_GayBerne = .false. |
| 723 |
|
FF_uses_EAM = .false. |
| 724 |
+ |
FF_uses_SC = .false. |
| 725 |
|
|
| 726 |
|
call getMatchingElementList(atypes, "is_Directional", .true., & |
| 727 |
|
nMatches, MatchList) |
| 738 |
|
call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList) |
| 739 |
|
if (nMatches .gt. 0) FF_uses_EAM = .true. |
| 740 |
|
|
| 741 |
+ |
call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList) |
| 742 |
+ |
if (nMatches .gt. 0) FF_uses_SC = .true. |
| 743 |
|
|
| 744 |
+ |
|
| 745 |
|
haveSaneForceField = .true. |
| 746 |
|
|
| 747 |
|
if (FF_uses_EAM) then |
| 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, & |
| 997 |
|
if (skipThisPair(atom1, atom2)) cycle inner |
| 998 |
|
|
| 999 |
|
if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then |
| 1000 |
< |
d_atm(1:3) = d_grp(1:3) |
| 1000 |
> |
d_atm(1) = d_grp(1) |
| 1001 |
> |
d_atm(2) = d_grp(2) |
| 1002 |
> |
d_atm(3) = d_grp(3) |
| 1003 |
|
ratmsq = rgrpsq |
| 1004 |
|
else |
| 1005 |
|
#ifdef IS_MPI |
| 1032 |
|
d_grp, rgrp, rCut) |
| 1033 |
|
#endif |
| 1034 |
|
vij = vij + vpair |
| 1035 |
< |
fij(1:3) = fij(1:3) + fpair(1:3) |
| 1035 |
> |
fij(1) = fij(1) + fpair(1) |
| 1036 |
> |
fij(2) = fij(2) + fpair(2) |
| 1037 |
> |
fij(3) = fij(3) + fpair(3) |
| 1038 |
|
endif |
| 1039 |
|
enddo inner |
| 1040 |
|
enddo |
| 1394 |
|
real( kind = dp ) :: d(3), scaled(3) |
| 1395 |
|
integer i |
| 1396 |
|
|
| 1397 |
< |
d(1:3) = q_j(1:3) - q_i(1:3) |
| 1397 |
> |
d(1) = q_j(1) - q_i(1) |
| 1398 |
> |
d(2) = q_j(2) - q_i(2) |
| 1399 |
> |
d(3) = q_j(3) - q_i(3) |
| 1400 |
|
|
| 1401 |
|
! Wrap back into periodic box if necessary |
| 1402 |
|
if ( SIM_uses_PBC ) then |
| 1419 |
|
else |
| 1420 |
|
! calc the scaled coordinates. |
| 1421 |
|
|
| 1405 |
– |
do i = 1, 3 |
| 1406 |
– |
scaled(i) = d(i) * HmatInv(i,i) |
| 1422 |
|
|
| 1423 |
< |
! wrap the scaled coordinates |
| 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) - dnint(scaled(1)) |
| 1430 |
> |
scaled(2) = scaled(2) - dnint(scaled(2)) |
| 1431 |
> |
scaled(3) = scaled(3) - dnint(scaled(3)) |
| 1432 |
|
|
| 1433 |
< |
scaled(i) = scaled(i) - anint(scaled(i)) |
| 1433 |
> |
! calc the wrapped real coordinates from the wrapped scaled |
| 1434 |
> |
! coordinates |
| 1435 |
|
|
| 1436 |
< |
! calc the wrapped real coordinates from the wrapped scaled |
| 1437 |
< |
! coordinates |
| 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 |
|
|
| 1415 |
– |
d(i) = scaled(i)*Hmat(i,i) |
| 1416 |
– |
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 |
|
|
