1C
2C CDDL HEADER START
3C
4C The contents of this file are subject to the terms of the Common Development
5C and Distribution License Version 1.0 (the "License").
6C
7C You can obtain a copy of the license at
8C http://www.opensource.org/licenses/CDDL-1.0.  See the License for the
9C specific language governing permissions and limitations under the License.
10C
11C When distributing Covered Code, include this CDDL HEADER in each file and
12C include the License file in a prominent location with the name LICENSE.CDDL.
13C If applicable, add the following below this CDDL HEADER, with the fields
14C enclosed by brackets "[]" replaced with your own identifying information:
15C
16C Portions Copyright (c) [yyyy] [name of copyright owner]. All rights reserved.
17C
18C CDDL HEADER END
19C
20
21C
22C Copyright (c) 2012, Regents of the University of Minnesota.  All rights reserved.
23C
24C Contributors:
25C    Valeriu Smirichinski
26C
27
28C*******************************************************************************
29C
30C  Release: This file is part of the openkim-api-v1.1.1 package.
31C
32C*******************************************************************************
33
34#include "KIM_API_status.h"
35
36C*******************************************************************************
37C
38C  FORTRAN 77 SUBROUTINE for calculating Lennard-Jones potential
39C
40C*******************************************************************************
41      SUBROUTINE ljpot_f77(sigma,epsilon,rr,v,dv)
42      IMPLICIT NONE
43
44C-------Transferred variables
45      DOUBLE PRECISION sigma,epsilon,rr,v,dv
46
47C-------Local variables
48      DOUBLE PRECISION sor2,sor6,sor12
49
50      IF (rr.LT.1.0E-14) STOP 'rr is zero'
51      sor2  = (sigma*sigma)/rr
52      sor6  = sor2*sor2*sor2
53      sor12 = sor6*sor6
54      v  = 4.D0*epsilon*(sor12 - sor6)
55      dv = 24.D0*epsilon*(-2.d0*sor12 + sor6)/sqrt(rr)
56      RETURN
57      END
58
59C*******************************************************************************
60C
61C  FORTRAN 77 SUBROUTINE to compute energy and forces on atoms from the
62C  positions.
63C
64C*******************************************************************************
65      SUBROUTINE calculate(cutoff,sigma,epsilon,pkim,x,f,ea,natom,
66     &                     ncontrib,en,e_flag,f_flag,eper_flag,
67     &                     kimget_neigh,ier)
68      IMPLICIT NONE
69
70C-------Transferred variables
71      DOUBLE PRECISION cutoff,sigma,epsilon
72      INTEGER kim
73      pointer (pkim,kim)
74      DOUBLE PRECISION x(3,natom),f(3,natom),ea(natom)
75      INTEGER natom,ier,e_flag,f_flag,eper_flag
76      DOUBLE PRECISION en
77      INTEGER kimget_neigh
78      EXTERNAL kimget_neigh
79
80C-------Local variables
81      DOUBLE PRECISION vij,dv,sumv,cut2,energycutoff
82      INTEGER i,j,jj,numnei,ncontrib
83      DOUBLE PRECISION r2
84      DOUBLE PRECISION r
85      DOUBLE PRECISION xi(3),xj(3),dx(3)
86      DOUBLE PRECISION Rij(3,512)
87      pointer (pRij,Rij)
88      INTEGER nei1atom(1)
89      pointer (pnei1atom,nei1atom)
90      INTEGER retcode,mode,request,atom
91
92      atom = 0
93      numnei = 0
94      IF (f_flag.EQ.1) THEN
95         DO 100 i = 1,natom
96            f(1,i) = 0.D0
97            f(2,i) = 0.D0
98            f(3,i) = 0.D0
99 100     CONTINUE
100      ENDIF
101      IF (eper_flag.EQ.1) THEN
102         DO 110 i = 1,natom
103 110        ea(i) = 0.D0
104      ENDIF
105      sumv = 0.D0
106      cut2 = cutoff*cutoff
107      CALL ljpot_f77(sigma,epsilon,cut2,energycutoff,dv)
108
109C     Iterator mode (reset iterator to beginning)
110      mode = 0
111      request = 0
112      retcode=kimget_neigh(pkim,mode,request,atom,
113     &                     numnei,pnei1atom,pRij)
114      IF (retcode.NE.KIM_STATUS_NEIGH_ITER_INIT_OK) THEN
115         WRITE(*,'("model_Ne_pure_LJ_NEIGH_PURE_*.F error: ",I5)')
116     &      retcode
117         ier = retcode
118         RETURN
119      ENDIF
120      retcode = KIM_STATUS_OK
121
122C     Loop over atoms and compute energy and forces
123 120  CONTINUE
124      IF (retcode.NE.KIM_STATUS_OK) GOTO 140
125
126C     Increment iterator
127      mode = 0
128      request = 1
129      retcode = kimget_neigh(pkim,mode,request,
130     &                       atom,numnei,pnei1atom,pRij)
131      IF (retcode.EQ.KIM_STATUS_NEIGH_ITER_PAST_END) THEN
132         ier = KIM_STATUS_OK
133         GOTO 140
134      ENDIF
135      IF (retcode.LT.KIM_STATUS_OK) THEN
136         WRITE(*,'("neigh iterator error:retcode",I5)'),retcode
137         ier = retcode
138         RETURN
139      ENDIF
140      i = atom
141      xi(1) = x(1,i)
142      xi(2) = x(2,i)
143      xi(3) = x(3,i)
144      DO 130 jj = 1,numnei
145         j = nei1atom(jj)
146         xj(1) = x(1,j)
147         xj(2) = x(2,j)
148         xj(3) = x(3,j)
149         dx(1) = xj(1)-xi(1)
150         dx(2) = xj(2)-xi(2)
151         dx(3) = xj(3)-xi(3)
152         r2 = dx(1)*dx(1) + dx(2)*dx(2) + dx(3)*dx(3)
153         IF (r2.LT.1.D-14)
154     &      WRITE(*,'("i=",I5,", j=",I5)') i,j
155         IF (r2.LE.cut2) THEN
156            CALL ljpot_f77(sigma,epsilon,r2,vij,dv)
157            IF (j.LE.ncontrib) THEN
158               sumv = sumv + vij-energycutoff
159            ELSE
160               sumv = sumv + 0.5d0*(vij-energycutoff)
161               dv = 0.5d0*dv;
162            ENDIF
163            IF (eper_flag.EQ.1) THEN
164               ea(i)=ea(i)+(vij-energycutoff)/2.D0
165               IF (j.LE.ncontrib) THEN
166                  ea(j)=ea(j)+(vij-energycutoff)/2.D0
167               ENDIF
168            ENDIF
169            IF (f_flag.EQ.1) THEN
170               r = sqrt(r2)
171               f(1,i) = f(1,i) + dv*dx(1)/r
172               f(2,i) = f(2,i) + dv*dx(2)/r
173               f(3,i) = f(3,i) + dv*dx(3)/r
174               f(1,j) = f(1,j) - dv*dx(1)/r
175               f(2,j) = f(2,j) - dv*dx(2)/r
176               f(3,j) = f(3,j) - dv*dx(3)/r
177            ENDIF
178         ENDIF
179 130  CONTINUE
180      GOTO 120
181 140  CONTINUE
182      IF (e_flag.EQ.1) en=sumv
183      RETURN
184      END
185