1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2! Copyright 2010. Los Alamos National Security, LLC. This material was ! 3! produced under U.S. Government contract DE-AC52-06NA25396 for Los Alamos ! 4! National Laboratory (LANL), which is operated by Los Alamos National ! 5! Security, LLC for the U.S. Department of Energy. The U.S. Government has ! 6! rights to use, reproduce, and distribute this software. NEITHER THE ! 7! GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC MAKES ANY WARRANTY, ! 8! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS ! 9! SOFTWARE. If software is modified to produce derivative works, such ! 10! modified software should be clearly marked, so as not to confuse it ! 11! with the version available from LANL. ! 12! ! 13! Additionally, this program is free software; you can redistribute it ! 14! and/or modify it under the terms of the GNU General Public License as ! 15! published by the Free Software Foundation; version 2.0 of the License. ! 16! Accordingly, this program is distributed in the hope that it will be ! 17! useful, but WITHOUT ANY WARRANTY; without even the implied warranty of ! 18! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General ! 19! Public License for more details. ! 20!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 21 22SUBROUTINE PAIRPOTNONEB 23 24 USE CONSTANTS_MOD 25 USE SETUPARRAY 26 USE PPOTARRAY 27 USE NEBLISTARRAY 28 USE VIRIALARRAY 29 USE MYPRECISION 30 31 IMPLICIT NONE 32 33 INTEGER :: I, NEWJ, J, K, PPSEL, BREAKLOOP 34 INTEGER :: PBCI, PBCJ, PBCK 35 REAL(LATTEPREC) :: JR1, JRCUT, R1, RCUT2, RCUT 36 REAL(LATTEPREC) :: FORCE, DC(3), RIJ(3) 37 REAL(LATTEPREC) :: MYR, MYR2, MYR3, MYR4, MAGR2, MAGR 38 REAL(LATTEPREC) :: UNIVPHI, JOINPHI, VDWPHI, CUTPHI, TMP 39 REAL(LATTEPREC) :: VIRUNIV(6), VIRJOIN(6), VIRVDW(6), VIRCUT(6) 40 REAL(LATTEPREC) :: FUNIV(3), FJOIN(3), FVDW(3), FCUT(3) 41 REAL(LATTEPREC) :: PHI, DPHI(3), EXPTMP, R6, FTMP(3) 42 REAL(LATTEPREC) :: POLYNOM, DPOLYNOM 43 IF (EXISTERROR) RETURN 44 45 ! 46 ! In this subroutine we add contributions in a strange way to ensure 47 ! numerical accuracy when switching between single and double precision. 48 ! If we don't do this, we get errors associated with adding very small 49 ! numbers to very large ones, and energies can be off by 0.01% or more. 50 ! 51 52 ! 53 ! There are 4 different parts to the pair potential: 54 ! 55 ! 1) Short range repulsion fitting to give bond lengths etc 56 ! 2) The joining function from JOINR1 TO JOINRCUT 57 ! 3) The vdW-type pair potential from JOINCUT to PPR1 58 ! 4) The final cut off tail from PPR1 TO PPRCUT 59 ! 60 61 UNIVPHI = ZERO 62 CUTPHI = ZERO 63 64 VIRUNIV = ZERO 65 VIRCUT = ZERO 66 67 IF (PPOTON .EQ. 1) THEN 68 69 DO I = 1, NATS 70 71 FUNIV = ZERO 72 FCUT = ZERO 73 74 ! Loop over all neighbors of I 75 76 DO J = 1, NATS 77 78 ! DO NEWJ = 1, TOTNEBPP(I) 79 80 ! J = NEBPP(1, NEWJ, I) 81 82 ! PBCI = NEBPP(2, NEWJ, I) 83 ! PBCJ = NEBPP(3, NEWJ, I) 84 ! PBCK = NEBPP(4, NEWJ, I) 85 86 DO K = 1, NOPPS 87 88 IF ((ATELE(I) .EQ. PPELE1(K) .AND. ATELE(J) .EQ. PPELE2(K)) & 89 .OR. (ATELE(J) .EQ. PPELE1(K) .AND. & 90 ATELE(I) .EQ. PPELE2(K))) THEN 91 92 PPSEL = K 93 94 R1 = POTCOEF(9,PPSEL) 95 RCUT = POTCOEF(10,PPSEL) 96 RCUT2 = RCUT*RCUT 97 98 ENDIF 99 100 ENDDO 101 102 RIJ(1) = CR(1,J) - CR(1,I) 103 RIJ(2) = CR(2,J) - CR(2,I) 104 RIJ(3) = CR(3,J) - CR(3,I) 105 106 ! RIJ(1) = CR(1,J) + REAL(PBCI)*BOX(1,1) + REAL(PBCJ)*BOX(2,1) + & 107 ! REAL(PBCK)*BOX(3,1) - CR(1,I) 108 109 ! RIJ(2) = CR(2,J) + REAL(PBCI)*BOX(1,2) + REAL(PBCJ)*BOX(2,2) + & 110 ! REAL(PBCK)*BOX(3,2) - CR(2,I) 111 112 ! RIJ(3) = CR(3,J) + REAL(PBCI)*BOX(1,3) + REAL(PBCJ)*BOX(2,3) + & 113 ! REAL(PBCK)*BOX(3,3) - CR(3,I) 114 115 MAGR2 = RIJ(1)*RIJ(1) + RIJ(2)*RIJ(2) + RIJ(3)*RIJ(3) 116 117 IF (MAGR2 .LE. RCUT2 .AND. J .NE. I) THEN 118 119 MAGR = SQRT(MAGR2) 120 121 ! Direction cosines 122 123 DC = RIJ/MAGR 124 125 IF (MAGR .LT. R1) THEN 126 127 ! CALL DUNIVSCALE(MAGR, POTCOEF(:,PPSEL), DC, PHI, DPHI) 128 129 POLYNOM = MAGR*(POTCOEF(2,PPSEL) + MAGR*(POTCOEF(3,PPSEL) + & 130 MAGR*(POTCOEF(4,PPSEL) + MAGR*POTCOEF(5,PPSEL)))) 131 132 PHI = POTCOEF(1,PPSEL)*EXP(POLYNOM) 133 134 DPOLYNOM = POTCOEF(2,PPSEL) + MAGR*(TWO*POTCOEF(3,PPSEL) + & 135 MAGR*(THREE*POTCOEF(4,PPSEL) + & 136 FOUR*POTCOEF(5,PPSEL)*MAGR)) 137 138 DPHI = -DC*PHI*DPOLYNOM 139 140 ! Hack! 141 EXPTMP = POTCOEF(6,PPSEL)*& 142 EXP( POTCOEF(7,PPSEL) * (MAGR - POTCOEF(8,PPSEL)) ) 143 ! R6 = MAGR2*MAGR2*MAGR2 144 145 ! UNIVPHI = UNIVPHI + PHI + EXPTMP - POTCOEF(8,PPSEL)/R6 146 147 UNIVPHI = UNIVPHI + PHI + EXPTMP 148 149 FTMP = DC*POTCOEF(7,PPSEL)*EXPTMP 150 ! FTMP = DC*(POTCOEF(7,PPSEL)*EXPTMP + & 151 ! SIX*POTCOEF(8,PPSEL)/(MAGR*R6)) 152 153 FUNIV = FUNIV - DPHI + FTMP 154 155 VIRUNIV(1) = VIRUNIV(1) - RIJ(1)*(DPHI(1) - FTMP(1)) 156 VIRUNIV(2) = VIRUNIV(2) - RIJ(2)*(DPHI(2) - FTMP(2)) 157 VIRUNIV(3) = VIRUNIV(3) - RIJ(3)*(DPHI(3) - FTMP(3)) 158 VIRUNIV(4) = VIRUNIV(4) - RIJ(1)*(DPHI(2) - FTMP(2)) 159 VIRUNIV(5) = VIRUNIV(5) - RIJ(2)*(DPHI(3) - FTMP(3)) 160 VIRUNIV(6) = VIRUNIV(6) - RIJ(3)*(DPHI(1) - FTMP(1)) 161 162 ELSE 163 164 MYR = MAGR - R1 165 166 CUTPHI = CUTPHI + POTCOEF(11,PPSEL) + & 167 MYR*(POTCOEF(12,PPSEL) + MYR*(POTCOEF(13,PPSEL) + & 168 MYR*(POTCOEF(14,PPSEL) + MYR*(POTCOEF(15,PPSEL) + & 169 MYR*POTCOEF(16,PPSEL))))) 170 171 FORCE = POTCOEF(12,PPSEL) + MYR*(TWO*POTCOEF(13,PPSEL) + & 172 MYR*(THREE*POTCOEF(14,PPSEL) + & 173 MYR*(FOUR*POTCOEF(15,PPSEL) + & 174 MYR*FIVE*POTCOEF(16,PPSEL)))) 175 176 FCUT = FCUT + DC*FORCE 177 178 VIRCUT(1) = VIRCUT(1) + RIJ(1)*DC(1)*FORCE 179 VIRCUT(2) = VIRCUT(2) + RIJ(2)*DC(2)*FORCE 180 VIRCUT(3) = VIRCUT(3) + RIJ(3)*DC(3)*FORCE 181 VIRCUT(4) = VIRCUT(4) + RIJ(1)*DC(2)*FORCE 182 VIRCUT(5) = VIRCUT(5) + RIJ(2)*DC(3)*FORCE 183 VIRCUT(6) = VIRCUT(6) + RIJ(3)*DC(1)*FORCE 184 185 ENDIF 186 187 ENDIF 188 189 ENDDO 190 191 FPP(1,I) = FUNIV(1) + FCUT(1) 192 FPP(2,I) = FUNIV(2) + FCUT(2) 193 FPP(3,I) = FUNIV(3) + FCUT(3) 194 195 ENDDO 196 197 EREP = HALF*(UNIVPHI + CUTPHI) 198 199 ! PRINT*, "EREP ", EREP 200 VIRPAIR = HALF*(VIRUNIV + VIRCUT) 201 202 ELSE 203 204 FPP = ZERO 205 EREP = ZERO 206 VIRPAIR = ZERO 207 208 ENDIF 209 210 RETURN 211 212END SUBROUTINE PAIRPOTNONEB 213 214