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