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