1 subroutine inisw 2 implicit real (a-h,o-z) 3 logical usecut,usesw 4 common /limits/ cutoff, cutof2,cuton, cuton2,conof3,usecut,usesw 5 6 cutof2 = cutoff*cutoff 7 cuton = cutoff - 1.0e0 8 9 cuton2 = cuton*cuton 10 tmp = cutof2 - cuton2 11 conof3 = tmp*tmp*tmp 12 isw = 1 13 14 15 return 16 end 17 18 real function swvdw(rij2) 19 implicit real (a-h,o-z) 20 logical usecut,usesw 21 common /limits/ cutoff, cutof2,cuton, cuton2,conof3,usecut,usesw 22 23 if (rij2.ge.cutof2) then 24 swvdw = 0.0e0 25 elseif (rij2.le.cuton2) then 26 swvdw = 1.0e0 27 else 28 rd1 = cutof2 - rij2 29 rd2 = cutof2 + 2.0e0*rij2 - 3.0e0*cuton2 30 swvdw = (rd1*rd1*rd2)/conof3 31 endif 32 33 return 34 end 35 36 real function swdvdw(rij2) 37c 38c derivative of VDW switch function 39c 40 implicit real (a-h,o-z) 41 logical usecut,usesw 42 common /limits/ cutoff, cutof2,cuton, cuton2,conof3,usecut,usesw 43 44 if (rij2.ge.cutof2) then 45 swdvdw = 0.0e0 46 elseif (rij2.le.cuton2) then 47 swdvdw = 1.0e0 48 else 49 rd1 = cutof2 - rij2 50c rd2 = rij2 - cuton2 51c swdvdw = 12.0e0*rd1*rd2/conof3 52 rd2 = cutof2 + 2.0e0*rij2 - 3.0e0*cuton2 53 rij = sqrt(rij2) 54 swdvdw = 4.0e0*rij*rd1*(rd1-rd2)/conof3 55 endif 56 57 return 58 end 59 60c swchg not used, instead same function as swvdw used for charge cutoff 61 62 subroutine swchg(rij2,s,ds) 63 implicit real (a-h,o-z) 64 logical usecut,usesw 65 common /limits/ cutoff, cutof2,cuton, cuton2,conof3,usecut,usesw 66 67 68 if (rij2.ge.cutof2) then 69 s = 0.0e0 70 ds = 0.0e0 71 else 72 frac = rij2/cutof2 73 frac2 = frac*frac 74 s = 1.0e0 - frac 75 s = s*s 76 ds = (frac2 - frac)*4.0e0 77c ds = ds/rij 78c we put the rij of de = de / r also in here 79c de = e*ds - e*s/r , de = de / r => 80c de = e*ds/rij - e*s/rij2 81c now de = e*ds - e*s 82 s = s/rij2 83 ds = ds/rij2 84 endif 85 86 return 87 end 88 89 subroutine bldlsd(coo,iresid,nlst,lst,istat) 90 implicit real (a-h,o-z), integer (i-n) 91 parameter (mxneib=200) 92 parameter (numres=50000) 93 logical usecut,usesw 94 common /limits/ cutoff, cutof2,cuton, cuton2,conof3,usecut,usesw 95 common /athlp/ iatoms, mxnat 96 common /residu/ qtot,ihsres, 97 & ires(numres),ibeg(numres),iend(numres) 98 logical box,cell,fast 99 common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast 100 logical doit 101 dimension coo(3,*),iresid(*),nlst(*),lst(*),vr(3),loclst(mxneib) 102 103 if (istat.eq.0) then 104 usesw = .false. 105 print*,'Not enough memory to use cutoff/switch' 106 return 107 endif 108 109 cutf1 = cutoff + 1.0e0 110 cutf2 = cutf1*cutf1 111 112c Build nonbonded interaction list, iatoms*mxneib 113c each atom can have a max of mxneib neighbours within cutoff of 8 angs 114c now not an atom is stored, but a residue 115 116 do i=1,iatoms 117 nlst(i) = 0 118 do j=i+1,iatoms 119 120 if (box) then 121 do k=1,3 122 vr(k) = coo(k,i) - coo(k,j) 123 end do 124 call reddis(vr) 125 r2 = vr(1)*vr(1) + vr(2)*vr(2) + vr(3)*vr(3) 126 doit = (r2.lt.cutf2) 127 else 128 doit = (dist2(coo(1,i),coo(1,j)).lt.cutf2) 129 endif 130 131 if (doit) then 132 if (nlst(i).lt.mxneib) then 133 ido = 1 134 do k=1,nlst(i) 135 if (iresid(j).eq.lst((i-1)*mxneib+k)) ido = 0 136 end do 137 if (ido.eq.1) then 138 nlst(i) = nlst(i) + 1 139 lst((i-1)*mxneib+nlst(i)) = iresid(j) 140 endif 141 else 142 print*,"neighbour list full for atom ",i 143 print*,"increase mxneib " 144 endif 145 endif 146 end do 147 end do 148 149 do i=1,ihsres 150 if (ires(i).gt.0) then 151 nloc = nlst(ibeg(i)) 152 do j=1,nloc 153 loclst(j) = lst((ibeg(i)-1)*mxneib+j) 154 end do 155 do j=ibeg(i)+1,iend(i) 156 do k=1,nlst(j) 157 ir = lst((j-1)*mxneib+k) 158 ido = 1 159 do l=1,nloc 160 if (ir.eq.loclst(l)) ido = 0 161 end do 162 if (ido.eq.1) then 163 if (nloc.lt.mxneib) then 164 nloc = nloc + 1 165 loclst(nloc) = ir 166 else 167 print*,"neighbour list full for atom ",ibeg(i) 168 print*,"increase mxneib " 169 endif 170 endif 171 end do 172 end do 173 do j=ibeg(i),iend(i) 174 nlst(j) = nloc 175 do l=1,nloc 176 lst((j-1)*mxneib+l) = loclst(l) 177 end do 178 end do 179 endif 180 end do 181 182 return 183 end 184 185 logical function resinc(iatom,nlst,lst,ires) 186 implicit real (a-h,o-z) 187 parameter (mxneib=200) 188 dimension lst(*) 189 190 191 resinc = .false. 192 do i=1,nlst 193 if (lst((iatom-1)*mxneib + i).eq.ires) then 194 resinc = .true. 195 return 196 endif 197 end do 198 199 return 200 end 201 202 real function dist2(a,b) 203c 204c determine distances between neighboring atoms 205c 206 implicit real (a-h,o-z) 207 dimension a(3) 208 dimension b(3) 209 210 d1 = a(1)-b(1) 211 d2 = a(2)-b(2) 212 d3 = a(3)-b(3) 213 214 dist2 = d1*d1 + d2*d2 + d3*d3 215 216 return 217 end 218 219 subroutine premul 220 implicit real (a-h,o-z), integer (i-n) 221 parameter (mxcls=50) 222 common /clspar/ ambvdwr(mxcls),ambvdwe(mxcls),mapagf(mxcls) 223 common /prem/ rs(mxcls,mxcls),es(mxcls,mxcls),iprem 224 225 do i=1,mxcls 226 do j=1,mxcls 227 rst = ambvdwr(i) + ambvdwr(j) 228 rs(i,j) = rst*rst*rst*rst*rst*rst 229 es(i,j) = sqrt(ambvdwe(i)*ambvdwe(j)) 230 end do 231 end do 232 233 call allerf 234 call alltab 235 236 iprem = 1 237 238 return 239 end 240 241 subroutine prmulr 242 implicit real(a-h,o-z), integer (i-n) 243 parameter (mxcls=50) 244 real ambvdwr,ambvdwe 245 common /clspar/ ambvdwr(mxcls),ambvdwe(mxcls),mapagf(mxcls) 246 common /premr/ rs(mxcls,mxcls),es(mxcls,mxcls),iprem 247 248 do i=1,mxcls 249 do j=1,mxcls 250 rst = real(ambvdwr(i) + ambvdwr(j)) 251 rs(i,j) = rst*rst*rst*rst*rst*rst 252 es(i,j) = sqrt(real(ambvdwe(i)*ambvdwe(j))) 253 end do 254 end do 255 256 iprem = 1 257 258 call alltab 259 260 return 261 end 262 263 subroutine alltab 264 implicit real (a-h,o-z), integer (i-n) 265 parameter (mxion=2000) 266 logical box,cell,fast 267 common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast 268 common /h2oer/ numwat,natnow,nion,iontyp,ionpl,iw(mxion),niwat 269 common /atab/ ialtab 270 271 if (ialtab.eq.1) return 272 273 call watcnt 274 call watlst(niwat) 275 call stutab 276 277 ialtab = 1 278 279 return 280 end 281 282 subroutine stutad(woo,whh,woh,dewoo,dewhh,dewoh) 283 implicit real (a-h,o-z), integer (i-n) 284 dimension woo(*),whh(*),woh(*),dewoo(*),dewhh(*),dewoh(*) 285 286 rcut = 9.0e0 287 rcinv = 1.0e0 / rcut 288 rc2inv = rcinv * rcinv 289 290 econv = 332.05382e0 291 qo = -0.8340e0 292 qh = 0.4170e0 293 vro = 1.7683e0 294 veo = 0.1520e0 295 vrh = 0.0001e0 296 veh = 0.0e0 297 298 poo = econv*qo*qo 299 rst = 2.0e0*vro 300 rst2 = rst*rst 301 pr6oo = rst2*rst2*rst2 302 epsoo = veo 303 304 phh = econv*qh*qh 305 pr6hh = 0.0e0 306 epshh = 0.0e0 307 308 poh = econv*qo*qh 309 rst = vro + vrh 310 rst2 = rst*rst 311 pr6oh = rst2*rst2*rst2 312 epsoh = sqrt(veo*veh) 313 314 nsize = 9000 315 316 do i=1,nsize 317 318 r = dble(i) / 1000.0e0 319 rinv = 1.0e0 / r 320 r2inv = rinv*rinv 321 r6inv = r2inv*r2inv*r2inv 322 323 p6 = pr6oo*r6inv 324 p12 = p6*p6 325 eq = poo * (rinv - rcinv + (r-rcut)*rc2inv) 326 ev = epsoo*(p12-2.0e0*p6) 327 328 woo(i) = eq + ev 329 330c if (i.lt.8000) print*,i,' woo ',woo(i), 331c & ' dwoo ',woo(i-1)-woo(i) 332 de = epsoo * (p12 - p6) * (-12.0e0) 333 dewoo(i) = (de-eq)*r2inv + poo*rinv*rc2inv 334 335 eq = phh * (rinv - rcinv + (r-rcut)*rc2inv) 336 337 whh(i) = eq 338 339 dewhh(i) = -eq*r2inv + phh*rinv*rc2inv 340 341 p6 = pr6oh*r6inv 342 p12 = p6*p6 343 eq = poh * (rinv - rcinv + (r-rcut)*rc2inv) 344 ev = epsoh*(p12-2.0e0*p6) 345 346 woh(i) = eq + ev 347 348 de = epsoh * (p12 - p6) * (-12.0e0) 349 dewoh(i) = (de-eq)*r2inv + poh*rinv*rc2inv 350 351 end do 352 353 do i=9000,10000 354 woo(i) = 0.0e0 355 woh(i) = 0.0e0 356 whh(i) = 0.0e0 357 dewoo(i) = 0.0e0 358 dewoh(i) = 0.0e0 359 dewhh(i) = 0.0e0 360 end do 361 362 return 363 end 364 365 subroutine etutad(nsize,woo,whh,woh,dewoo,dewhh,dewoh) 366 implicit real (a-h,o-z), integer (i-n) 367 dimension woo(*),whh(*),woh(*),dewoo(*),dewhh(*),dewoh(*) 368 369 econv = 332.05382e0 370 qo = -0.8340e0 371 qh = 0.4170e0 372 vro = 1.7683e0 373 veo = 0.1520e0 374 vrh = 0.0001e0 375 veh = 0.0e0 376 377 poo = econv*qo*qo 378 rst = 2.0e0*vro 379 rst2 = rst*rst 380 pr6oo = rst2*rst2*rst2 381 epsoo = veo 382 383 phh = econv*qh*qh 384 pr6hh = 0.0e0 385 epshh = 0.0e0 386 387 poh = econv*qo*qh 388 rst = vro + vrh 389 rst2 = rst*rst 390 pr6oh = rst2*rst2*rst2 391 epsoh = sqrt(veo*veh) 392 393 do i=1,nsize 394 395 r = real(i) / 1000.0e0 396 rinv = 1.0e0 / r 397 r2inv = rinv*rinv 398 r6inv = r2inv*r2inv*r2inv 399 400 p6 = pr6oo*r6inv 401 p12 = p6*p6 402 eq = poo * rinv 403 ev = epsoo*(p12-2.0e0*p6) 404 405 woo(i) = eq + ev 406 407 de = epsoo * (p12 - p6) * (-12.0e0) 408 dewoo(i) = (de-eq)*r2inv 409 410 eq = phh * rinv 411 412 whh(i) = eq 413 414 dewhh(i) = -eq*r2inv 415 416 p6 = pr6oh*r6inv 417 p12 = p6*p6 418 eq = poh * rinv 419 ev = epsoh*(p12-2.0e0*p6) 420 421 woh(i) = eq + ev 422 423 de = epsoh * (p12 - p6) * (-12.0e0) 424 dewoh(i) = (de-eq)*r2inv 425 426 end do 427 428 return 429 end 430 431 subroutine allerd(aprerf,aprexp) 432 implicit real (a-h,o-z), integer (i-n) 433 real aprerf,aprexp,alaf,pi,spi,f 434 dimension aprerf(*),aprexp(*) 435 436 alfa = 0.3e0 437 pi = 4.e0*atan(1.e0) 438 spi = sqrt(pi) 439 440 do i=1,1800000 441 442 r = dble(i) / 1000000.0e0 443 aprerf(i) = real(erfc(r)) 444 445 end do 446 447 do i=1,3240000 448 449 f = - float(i) / 1000000.0e0 450 aprexp(i) = 2.0*alfa*exp(f)/spi 451 452 end do 453 454 return 455 end 456 457c subroutine allerd(aprerf,aprexp) 458c implicit real (a-h,o-z), integer (i-n) 459c real aprerf,aprexp,alaf,pi,spi,f 460c dimension aprerf(*),aprexp(*) 461c 462c alfa = 0.2e0 463c pi = 4.e0*atan(1.e0) 464c spi = sqrt(pi) 465c 466c do i=1,3000000 467c 468c r = dble(i) / 1000000.0e0 469c aprerf(i) = real(erfc(r)) 470c 471c end do 472c 473c do i=1,9000000 474c 475c f = - float(i) / 1000000.0e0 476c aprexp(i) = 2.0*alfa*exp(f)/spi 477c 478c end do 479c 480c return 481c end 482 483 subroutine apperfd(r,fc,aprerf) 484 implicit real (a-h,o-z), integer (i-n) 485 real aprerf,r1000,x,w 486 dimension aprerf(*) 487 488 thou = 1000000.0e0 489 490 if (r.le.10.0e0) then 491 r1000 = real(r)*thou 492 ir = int(r1000) 493 x = (r1000 - float(ir)) 494 fc = dble((1.0e0-x)*aprerf(ir) + x*aprerf(ir+1)) 495 endif 496 497 return 498 end 499 500