1 logical function argos_prep_misfit(lfnout,lfnmat, 2 + isfnd,xs,csa,qsa,isgm,msa,nsa, 3 + idsb,cdsb,msb,nsb,grid,mgrid,ngrid,gdist,nmis,iwater,iwfnd,xw, 4 + qwa,mwm,mwa,nwm,nwa,npbtyp,box,fcount,nrgrid,iogrid,rogrid) 5c 6c $Id$ 7c 8 implicit none 9c 10#include "util.fh" 11c 12 logical md_zmat 13 external md_zmat 14c 15 integer lfnout,lfnmat 16 integer msa,nsa,msb,nsb,mgrid,nmis,iwater,mwm,mwa,nwm,nwa 17 integer isfnd(msa),idsb(2,msb),iwfnd(mwm),npbtyp,isgm(msa) 18 real*8 cdsb(msb),xs(3,msa),qsa(msa) 19 character*16 csa(msa) 20 real*8 xw(3,mwa,mwm),qwa(mwa),box(3),fcount 21 character*255 filmat 22 integer nrgrid,iogrid(5),rogrid(2,5) 23c 24 real*8 grid(4,mgrid),distx,disty,distz 25 integer i,ic,j,nn,nndx(10),nnn,nnndx(10),nm,miss 26 real*8 d(3),c(3),p(3),t(3),r,dr,dd,px(3),pn(3),pz(3),pw(3) 27 real*8 xmax(3),xmin(3),dx(3),dist2,gdist,gdist2,angle,boxh(3) 28 integer ix,iy,iz,k,l,m,imax,imin,ngrid,iwm,isa,iwm2,iwa,num 29c 30 integer isg,ioff,jz,kz,lz 31 real*8 dz,az,tz,dgrid 32 logical lmat,fndone 33c 34 if(nmis.gt.0) then 35c 36c process z-matrices for missing atoms 37c ------------------------------------ 38c 39 32 continue 40 dr=0.0d0 41 ioff=0 42 isg=0 43 lmat=.false. 44 fndone=.false. 45 do 30 i=1,nsa 46 if(isgm(i).ne.isg) then 47 isg=isgm(i) 48 ioff=i-1 49 filmat=csa(i)(1:index(csa(i),' ')-1)//'.mat ' 50 lmat=.false. 51 if(lmat) close(unit=lfnmat) 52 open(unit=lfnmat,file=filmat(1:index(filmat,' ')-1), 53 + form='formatted',status='old',err=30) 54 lmat=.true. 55 endif 56 if(lmat.and.isfnd(i).ne.1) then 57 rewind(lfnmat) 58 31 continue 59 read(lfnmat,1000,end=30,err=30) iz,jz,kz,lz,dz,az,tz 60 1000 format(4i5,3f12.6) 61 if(i.eq.iz+ioff.and.isfnd(ioff+jz).ne.0.and. 62 + isfnd(ioff+kz).ne.0.and.isfnd(ioff+lz).ne.0) then 63 if(md_zmat(xs(1,i),xs(1,ioff+jz),xs(1,ioff+kz),xs(1,ioff+lz), 64 + dz,az,tz)) then 65 isfnd(i)=1 66 write(*,'(a,i7,a,a,3f10.5)') ' Generated from z-matrix: ', 67 + isg,':',csa(i),(xs(k,i),k=1,3) 68 fndone=.true. 69 endif 70 else 71 goto 31 72 endif 73 endif 74 30 continue 75 if(lmat) close(unit=lfnmat) 76 if(fndone) goto 32 77c 78 if(util_print('restart',print_default)) then 79 write(lfnout,2000) 80 2000 format(' Z-matrix definitions done',/) 81 endif 82c 83c generate missing coordinates based on simple rules 84c -------------------------------------------------- 85c 86 gdist2=gdist*gdist 87 boxh(1)=0.5d0*box(1) 88 boxh(2)=0.5d0*box(2) 89 boxh(3)=0.5d0*box(3) 90c 91c loop over all solute atoms 92c -------------------------- 93c 94 ngrid=0 95 miss=0 96c 97 do 1 i=1,nsa 98c 99c for each atom without coordinates 100c --------------------------------- 101c 102 if(isfnd(i).ne.1) then 103 if(util_print('coordinates',print_debug)) then 104 write(lfnout,'(a,i7,5x,a16)') 'Not found: atom ',i,csa(i) 105 endif 106 nn=0 107c 108c find all bonded atoms 109c --------------------- 110c 111 do 2 j=1,nsb 112 if(idsb(1,j).eq.i) then 113 nn=nn+1 114 nndx(nn)=idsb(2,j) 115 endif 116 if(idsb(2,j).eq.i) then 117 nn=nn+1 118 nndx(nn)=idsb(1,j) 119 endif 120 2 continue 121c 122 if(nn.eq.0) goto 1 123c 124c if bound to single other existing atom 125c 126 if(nn.eq.1.and.isfnd(nndx(1)).eq.1) then 127c 128c find atoms bound to that atom 129c 130 nnn=0 131 do 3 j=1,nsb 132 if(idsb(1,j).eq.nndx(1)) then 133 nnn=nnn+1 134 nnndx(nnn)=idsb(2,j) 135 endif 136 if(idsb(2,j).eq.nndx(1)) then 137 nnn=nnn+1 138 nnndx(nnn)=idsb(1,j) 139 endif 140 if(idsb(1,j).eq.i.or.idsb(2,j).eq.i) dr=cdsb(j) 141 3 continue 142c 143c count the number of missing atoms 144c 145 nm=0 146 do 4 j=1,nnn 147 if(isfnd(nnndx(j)).ne.1) nm=nm+1 148 4 continue 149c 150c if j is only missing atom 151c 152 if(nm.eq.1) then 153 c(1)=xs(1,nndx(1)) 154 c(2)=xs(2,nndx(1)) 155 c(3)=xs(3,nndx(1)) 156 d(1)=0.0d0 157 d(2)=0.0d0 158 d(3)=0.0d0 159 num=0 160 do 5 j=1,nnn 161 if(nnndx(j).ne.i) then 162 d(1)=d(1)+xs(1,nnndx(j)) 163 d(2)=d(2)+xs(2,nnndx(j)) 164 d(3)=d(3)+xs(3,nnndx(j)) 165 num=num+1 166 endif 167 5 continue 168 d(1)=d(1)/dble(num)-c(1) 169 d(2)=d(2)/dble(num)-c(2) 170 d(3)=d(3)/dble(num)-c(3) 171 r=sqrt(d(1)*d(1)+d(2)*d(2)+d(3)*d(3)) 172 if(r.lt.0.0001) 173 + call md_abort('argos_prep_misfit: vector too small',9001) 174 xs(1,i)=c(1)-dr*d(1)/r 175 xs(2,i)=c(2)-dr*d(2)/r 176 xs(3,i)=c(3)-dr*d(3)/r 177 if(nnn.eq.2) then 178 p(1)=xs(1,i) 179 p(2)=xs(2,i) 180 p(3)=xs(3,i) 181 t(1)=p(2)*c(3)-c(2)*p(3)+c(1) 182 t(2)=p(3)*c(1)-c(3)*p(1)+c(2) 183 t(3)=p(1)*c(2)-c(1)*p(2)+c(3) 184 angle=1.047198d0 185 call rotate(c,t,angle,p,xs(1,i)) 186 endif 187 isfnd(i)=1 188 endif 189c 190c if j is one of two missing atoms 191c 192 if(nm.eq.2) then 193 c(1)=xs(1,nndx(1)) 194 c(2)=xs(2,nndx(1)) 195 c(3)=xs(3,nndx(1)) 196 d(1)=0.0d0 197 d(2)=0.0d0 198 d(3)=0.0d0 199 ic=0 200 num=0 201 do 6 j=1,nnn 202 if(nnndx(j).ne.i.and.isfnd(nnndx(j)).eq.1) then 203 d(1)=d(1)+xs(1,nnndx(j)) 204 d(2)=d(2)+xs(2,nnndx(j)) 205 d(3)=d(3)+xs(3,nnndx(j)) 206 ic=nnndx(j) 207 num=num+1 208 endif 209 6 continue 210 d(1)=d(1)/dble(num)-c(1) 211 d(2)=d(2)/dble(num)-c(2) 212 d(3)=d(3)/dble(num)-c(3) 213 r=sqrt(d(1)*d(1)+d(2)*d(2)+d(3)*d(3)) 214 if(r.lt.0.0001) 215 + call md_abort('argos_prep_misfit: vector too small',9002) 216 d(1)=dr*d(1)/r 217 d(2)=dr*d(2)/r 218 d(3)=dr*d(3)/r 219 if(nnn.eq.3) then 220 if(abs(d(1)).le.abs(d(2)).and.abs(d(1)).le.abs(d(3))) then 221 t(1)=d(1)+c(1) 222 t(2)=d(3)+c(2) 223 t(3)=-d(2)+c(3) 224 elseif(abs(d(2)).le.abs(d(1)).and.abs(d(2)).le.abs(d(3))) then 225 t(1)=d(3)+c(1) 226 t(2)=d(2)+c(2) 227 t(3)=-d(1)+c(3) 228 else 229 t(1)=d(2)+c(1) 230 t(2)=-d(1)+c(2) 231 t(3)=d(3)+c(3) 232 endif 233 d(1)=d(1)+c(1) 234 d(2)=d(2)+c(2) 235 d(3)=d(3)+c(3) 236 angle=2.094395102d0 237 call rotate(c,t,angle,d,p) 238 xs(1,i)=p(1) 239 xs(2,i)=p(2) 240 xs(3,i)=p(3) 241 isfnd(i)=1 242 else 243 d(1)=d(1)+c(1) 244 d(2)=d(2)+c(2) 245 d(3)=d(3)+c(3) 246 r=sqrt((xs(1,ic)-c(1))**2+(xs(2,ic)-c(2))**2+(xs(3,ic)-c(3))**2) 247 if(r.lt.0.0001) 248 + call md_abort('argos_prep_misfit: vector too small',9004) 249 t(1)=dr*(xs(1,ic)-c(1))/r+c(1) 250 t(2)=dr*(xs(2,ic)-c(2))/r+c(2) 251 t(3)=dr*(xs(3,ic)-c(3))/r+c(3) 252 angle=1.570796327d0 253 call rotate(c,d,angle,t,p) 254 d(1)=d(1)-c(1) 255 d(2)=d(2)-c(2) 256 d(3)=d(3)-c(3) 257 t(1)=t(1)-c(1) 258 t(2)=t(2)-c(2) 259 t(3)=t(3)-c(3) 260 c(1)=d(2)*t(3)-t(2)*d(3)+xs(1,nndx(1)) 261 c(2)=d(3)*t(1)-t(3)*d(1)+xs(2,nndx(1)) 262 c(3)=d(1)*t(2)-t(1)*d(2)+xs(3,nndx(1)) 263 d(1)=xs(1,nndx(1)) 264 d(2)=xs(2,nndx(1)) 265 d(3)=xs(3,nndx(1)) 266 angle=3.141592654d0 267 call rotate(d,c,angle,p,t) 268 xs(1,i)=t(1) 269 xs(2,i)=t(2) 270 xs(3,i)=t(3) 271 isfnd(i)=1 272 endif 273 endif 274c 275c if j is one of three missing atoms 276c 277 if(nm.eq.3) then 278 c(1)=xs(1,nndx(1)) 279 c(2)=xs(2,nndx(1)) 280 c(3)=xs(3,nndx(1)) 281 d(1)=0.0d0 282 d(2)=0.0d0 283 d(3)=0.0d0 284 ic=0 285 num=0 286 do 7 j=1,nnn 287 if(nnndx(j).ne.i.and.isfnd(nnndx(j)).eq.1) then 288 d(1)=d(1)+xs(1,nnndx(j)) 289 d(2)=d(2)+xs(2,nnndx(j)) 290 d(3)=d(3)+xs(3,nnndx(j)) 291 ic=nnndx(j) 292 num=num+1 293 endif 294 7 continue 295 d(1)=d(1)/dble(num)-c(1) 296 d(2)=d(2)/dble(num)-c(2) 297 d(3)=d(3)/dble(num)-c(3) 298 r=sqrt(d(1)*d(1)+d(2)*d(2)+d(3)*d(3)) 299 if(r.lt.0.0001) 300 + call md_abort('argos_prep_misfit: vector too small',9005) 301 d(1)=dr*d(1)/r 302 d(2)=dr*d(2)/r 303 d(3)=dr*d(3)/r 304 if(nnn.eq.4) then 305 if(abs(d(1)).le.abs(d(2)).and.abs(d(1)).le.abs(d(3))) then 306 t(1)=d(1)+c(1) 307 t(2)=d(3)+c(2) 308 t(3)=-d(2)+c(3) 309 elseif(abs(d(2)).le.abs(d(1)).and.abs(d(2)).le.abs(d(3))) then 310 t(1)=d(3)+c(1) 311 t(2)=d(2)+c(2) 312 t(3)=-d(1)+c(3) 313 else 314 t(1)=d(2)+c(1) 315 t(2)=-d(1)+c(2) 316 t(3)=d(3)+c(3) 317 endif 318 d(1)=d(1)+c(1) 319 d(2)=d(2)+c(2) 320 d(3)=d(3)+c(3) 321 angle=1.910633236d0 322 call rotate(c,t,angle,d,p) 323 xs(1,i)=p(1) 324 xs(2,i)=p(2) 325 xs(3,i)=p(3) 326 isfnd(i)=1 327 else 328c 329c three missing hydrogens from atom with other than 4 neighbors 330c 331 call md_abort('argos_prep_misfit: not implemented for atom',i) 332 endif 333 endif 334c 335 endif 336 endif 337c 338 1 continue 339c 340 if(util_print('restart',print_default)) then 341 write(lfnout,2001) 342 2001 format(' Geometric rules definitions done',/) 343 endif 344c 345 do 8 i=1,nsa 346c 347c for each atom without coordinates 348c --------------------------------- 349c 350 if(isfnd(i).ne.1) then 351 nn=0 352c 353c find all bonded atoms 354c --------------------- 355c 356 do 9 j=1,nsb 357 if(idsb(1,j).eq.i) then 358 nn=nn+1 359 nndx(nn)=idsb(2,j) 360 endif 361 if(idsb(2,j).eq.i) then 362 nn=nn+1 363 nndx(nn)=idsb(1,j) 364 endif 365 9 continue 366c 367c 368c if not bound to any atom 369c ------------------------ 370c 371 if(nn.eq.0) then 372 if(util_print('coordinates',print_debug)) then 373 write(lfnout,'(a,i7)') 'Setting up grid for atom ',i 374 endif 375c 376c setup grid 377c 378 if(ngrid.eq.0) then 379 do 10 j=1,nsa 380 if(isfnd(j).eq.1) then 381 xmax(1)=xs(1,j) 382 xmax(2)=xs(2,j) 383 xmax(3)=xs(3,j) 384 xmin(1)=xs(1,j) 385 xmin(2)=xs(2,j) 386 xmin(3)=xs(3,j) 387 goto 11 388 endif 389 10 continue 390 11 continue 391 do 12 k=1,3 392 do 13 j=1,nsa 393 if(isfnd(j).eq.1) then 394 if(xmax(k).lt.xs(k,j)) xmax(k)=xs(k,j) 395 if(xmin(k).gt.xs(k,j)) xmin(k)=xs(k,j) 396 endif 397 13 continue 398 dx(k)=(xmax(k)-xmin(k)+4.0d0*gdist)/(mgrid-1) 399 xmin(k)=xmin(k)-2.0d0*gdist 400 12 continue 401c 402 if(ngrid.gt.0.and.util_print('coordinates',print_debug)) then 403 do 555 k=1,nrgrid 404 write(lfnout,'(a,i5,2f12.5)') 'grid option ', 405 + iogrid(k),rogrid(1,k),rogrid(2,k) 406 555 continue 407 endif 408 ngrid=1 409 do 14 ix=1,mgrid 410 do 15 iy=1,mgrid 411 do 16 iz=1,mgrid 412 grid(1,ngrid)=xmin(1)+dble(ix-1)*dx(1) 413 grid(2,ngrid)=xmin(2)+dble(iy-1)*dx(2) 414 grid(3,ngrid)=xmin(3)+dble(iz-1)*dx(3) 415 if(nrgrid.gt.0) then 416 do 217 k=1,nrgrid 417 if(iogrid(k).eq.1) then 418 if(grid(1,ngrid).lt.rogrid(1,k)) goto 16 419 if(grid(1,ngrid).gt.rogrid(2,k)) goto 16 420 elseif(iogrid(k).eq.2) then 421 if(grid(1,ngrid).gt.rogrid(1,k).and.grid(1,ngrid).lt.rogrid(2,k)) 422 + goto 16 423 elseif(iogrid(k).eq.3) then 424 if(grid(2,ngrid).lt.rogrid(1,k)) goto 16 425 if(grid(2,ngrid).gt.rogrid(2,k)) goto 16 426 elseif(iogrid(k).eq.4) then 427 if(grid(2,ngrid).gt.rogrid(1,k).and.grid(2,ngrid).lt.rogrid(2,k)) 428 + goto 16 429 elseif(iogrid(k).eq.5) then 430 if(grid(3,ngrid).lt.rogrid(1,k)) goto 16 431 if(grid(3,ngrid).gt.rogrid(2,k)) goto 16 432 elseif(iogrid(k).eq.6) then 433 if(grid(3,ngrid).gt.rogrid(1,k).and.grid(3,ngrid).lt.rogrid(2,k)) 434 + goto 16 435 elseif(iogrid(k).eq.7) then 436 dgrid=sqrt(grid(1,ngrid)**2+grid(2,ngrid)**2) 437 if(dgrid.lt.rogrid(1,k)) goto 16 438 if(dgrid.gt.rogrid(2,k)) goto 16 439 elseif(iogrid(k).eq.8) then 440 dgrid=sqrt(grid(1,ngrid)**2+grid(2,ngrid)**2) 441 if(dgrid.gt.rogrid(1,k).and.dgrid.lt.rogrid(2,k)) goto 16 442 elseif(iogrid(k).eq.9) then 443 dgrid=sqrt(grid(1,ngrid)**2+grid(2,ngrid)**2+grid(3,ngrid)**2) 444 if(dgrid.lt.rogrid(1,k)) goto 16 445 if(dgrid.gt.rogrid(2,k)) goto 16 446 elseif(iogrid(k).eq.10) then 447 dgrid=sqrt(grid(1,ngrid)**2+grid(2,ngrid)**2+grid(2,ngrid)**2) 448 if(dgrid.gt.rogrid(1,k).and.dgrid.lt.rogrid(2,k)) goto 16 449 endif 450 217 continue 451 endif 452c write(*,'(a,4i5)') 'grid ',ix,iy,iz 453 do 17 k=1,nsa 454 if(isfnd(k).eq.1) then 455 if(npbtyp.eq.0) then 456 dist2=(grid(1,ngrid)-xs(1,k))**2+(grid(2,ngrid)-xs(2,k))**2+ 457 + (grid(3,ngrid)-xs(3,k))**2 458 else 459 distx=grid(1,ngrid)-xs(1,k) 460 disty=grid(2,ngrid)-xs(2,k) 461 distz=grid(3,ngrid)-xs(3,k) 462 if(abs(distx).gt.boxh(1)) distx=distx-nint(distx/box(1))*box(1) 463 if(abs(disty).gt.boxh(2)) disty=disty-nint(distx/box(2))*box(2) 464 if(abs(distz).gt.boxh(3)) distz=distz-nint(distx/box(3))*box(3) 465 dist2=distx**2+disty**2+distz**2 466 endif 467 if(dist2.lt.gdist2) goto 16 468 endif 469 17 continue 470 if(iwater.gt.0) then 471 do 117 k=1,nwm 472 if(iwfnd(k).gt.0) then 473 do 118 l=1,iwfnd(k) 474 if(npbtyp.eq.0) then 475 dist2=(grid(1,ngrid)-xw(1,l,k))**2+(grid(2,ngrid)-xw(2,l,k))**2+ 476 + (grid(3,ngrid)-xw(3,l,k))**2 477 else 478 distx=grid(1,ngrid)-xw(1,l,k) 479 disty=grid(2,ngrid)-xw(2,l,k) 480 distz=grid(3,ngrid)-xw(3,l,k) 481 if(abs(distx).gt.boxh(1)) distx=distx-nint(distx/box(1))*box(1) 482 if(abs(disty).gt.boxh(2)) disty=disty-nint(distx/box(2))*box(2) 483 if(abs(distz).gt.boxh(3)) distz=distz-nint(distx/box(3))*box(3) 484 dist2=distx**2+disty**2+distz**2 485 endif 486 if(dist2.lt.gdist2) goto 16 487 118 continue 488 endif 489 117 continue 490 endif 491 ngrid=ngrid+1 492 16 continue 493 15 continue 494 14 continue 495 ngrid=ngrid-1 496c 497 if(ngrid.le.0) call md_abort('No grid points 0',9999) 498c 499c calculate grid potential 500c 501 do 18 k=1,ngrid 502 grid(4,k)=0.0d0 503 do 19 j=1,nsa 504 if(isfnd(j).eq.1) then 505 if(npbtyp.eq.0) then 506 dist2=dsqrt((grid(1,k)-xs(1,j))**2+ 507 + (grid(2,k)-xs(2,j))**2+(grid(3,k)-xs(3,j))**2) 508 else 509 distx=grid(1,ngrid)-xs(1,k) 510 disty=grid(2,ngrid)-xs(2,k) 511 distz=grid(3,ngrid)-xs(3,k) 512 if(abs(distx).gt.boxh(1)) distx=distx-nint(distx/box(1))*box(1) 513 if(abs(disty).gt.boxh(2)) disty=disty-nint(distx/box(2))*box(2) 514 if(abs(distz).gt.boxh(3)) distz=distz-nint(distx/box(3))*box(3) 515 dist2=distx**2+disty**2+distz**2 516 endif 517 grid(4,k)=grid(4,k)+fcount*qsa(j)/dist2 518 endif 519 19 continue 520 if(iwater.gt.0) then 521 do 119 m=1,nwm 522 if(iwfnd(m).gt.0) then 523 do 120 l=1,iwfnd(m) 524 if(npbtyp.eq.0) then 525 dist2=(grid(1,ngrid)-xw(1,l,m))**2+(grid(2,ngrid)-xw(2,l,m))**2+ 526 + (grid(3,ngrid)-xw(3,l,m))**2 527 else 528 distx=grid(1,ngrid)-xw(1,l,m) 529 disty=grid(2,ngrid)-xw(2,l,m) 530 distz=grid(3,ngrid)-xw(3,l,m) 531 if(abs(distx).gt.boxh(1)) distx=distx-nint(distx/box(1))*box(1) 532 if(abs(disty).gt.boxh(2)) disty=disty-nint(distx/box(2))*box(2) 533 if(abs(distz).gt.boxh(3)) distz=distz-nint(distx/box(3))*box(3) 534 dist2=distx**2+disty**2+distz**2 535 endif 536 grid(4,k)=grid(4,k)+qwa(l)/dist2 537 120 continue 538 endif 539 119 continue 540 endif 541 18 continue 542c 543 endif 544c 545 if(ngrid.le.0) call md_abort('No grid points 1',9999) 546c 547c find extrema 548c 549 imax=1 550 imin=1 551 do 20 j=1,ngrid 552 if(grid(4,j).gt.grid(4,imax)) imax=j 553 if(grid(4,j).lt.grid(4,imin)) imin=j 554 20 continue 555c 556c choose best grid point 557c 558 if(qsa(i).gt.0.0d0) then 559 xs(1,i)=grid(1,imin) 560 xs(2,i)=grid(2,imin) 561 xs(3,i)=grid(3,imin) 562 else 563 xs(1,i)=grid(1,imax) 564 xs(2,i)=grid(2,imax) 565 xs(3,i)=grid(3,imax) 566 imin=imax 567 endif 568 isfnd(i)=1 569 if(util_print('restart',print_medium)) then 570 write(lfnout,6000) fcount*qsa(i)*grid(4,imin)*138.9354d0 571 6000 format(' Added counterion with energy ',t40,f12.3,' kJ/mol') 572 endif 573c 574c remove all grid points within gdist from imin from grid 575c 576 d(1)=grid(1,imin) 577 d(2)=grid(2,imin) 578 d(3)=grid(3,imin) 579c 580 imax=ngrid 581 ngrid=0 582 do 21 j=1,imax 583 if(npbtyp.eq.0) then 584 dist2=(grid(1,j)-d(1))**2+(grid(2,j)-d(2))**2+(grid(3,j)-d(3))**2 585 else 586 distx=grid(1,j)-d(1) 587 disty=grid(2,j)-d(2) 588 distz=grid(3,j)-d(3) 589 if(abs(distx).gt.boxh(1)) distx=distx-nint(distx/box(1))*box(1) 590 if(abs(disty).gt.boxh(2)) disty=disty-nint(distx/box(2))*box(2) 591 if(abs(distz).gt.boxh(3)) distz=distz-nint(distx/box(3))*box(3) 592 dist2=distx**2+disty**2+distz**2 593 endif 594 if(dist2.ge.gdist2) then 595 ngrid=ngrid+1 596 grid(1,ngrid)=grid(1,j) 597 grid(2,ngrid)=grid(2,j) 598 grid(3,ngrid)=grid(3,j) 599 dist2=dsqrt((grid(1,ngrid)-xs(1,i))**2+ 600 + (grid(2,ngrid)-xs(2,i))**2+(grid(3,ngrid)-xs(3,i))**2) 601 grid(4,ngrid)=grid(4,j)+fcount*qsa(i)/dist2 602 endif 603 21 continue 604c 605 if(ngrid.le.0) call md_abort('No grid points left',9999) 606c 607 else 608 miss=miss+1 609 if(util_print('restart',print_none)) then 610 write(lfnout,9998) isgm(i),csa(i) 611 9998 format(' Coordinates missing for atom ',i5,':',a) 612 endif 613 endif 614 endif 615 8 continue 616c 617 if(util_print('restart',print_default)) then 618 write(lfnout,2002) 619 2002 format(' Grid specifications for non-bonded atoms done',/) 620 endif 621c 622 if(miss.ne.0) then 623 call md_abort('argos_prep_misfit: missing non hydrogen',9999) 624 endif 625c 626 endif 627c 628c crystal waters 629c 630 if(iwater.gt.0) then 631 do 23 iwm=1,nwm 632 if(iwfnd(iwm).eq.1) then 633 px(1)=0.0d0 634 px(2)=0.0d0 635 px(3)=0.0d0 636 do 24 isa=1,nsa 637 dx(1)=(xs(1,isa)-xw(1,1,iwm)) 638 dx(2)=(xs(2,isa)-xw(2,1,iwm)) 639 dx(3)=(xs(3,isa)-xw(3,1,iwm)) 640 dd=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3) 641 px(1)=px(1)+qsa(isa)*dx(1)/dd 642 px(2)=px(2)+qsa(isa)*dx(2)/dd 643 px(3)=px(3)+qsa(isa)*dx(3)/dd 644 24 continue 645 do 25 iwm2=1,nwm 646 if(iwfnd(iwm2).eq.3) then 647 do 26 iwa=1,3 648 dx(1)=(xw(1,iwa,iwm2)-xw(1,1,iwm)) 649 dx(2)=(xw(2,iwa,iwm2)-xw(2,1,iwm)) 650 dx(3)=(xw(3,iwa,iwm2)-xw(3,1,iwm)) 651 dd=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3) 652 px(1)=px(1)+qwa(iwa)*dx(1)/dd 653 px(2)=px(2)+qwa(iwa)*dx(2)/dd 654 px(3)=px(3)+qwa(iwa)*dx(3)/dd 655 26 continue 656 endif 657 25 continue 658 dd=sqrt(px(1)*px(1)+px(2)*px(2)+px(3)*px(3)) 659 pw(1)=-0.1d0*px(1)/dd 660 pw(2)=-0.1d0*px(2)/dd 661 pw(3)=-0.1d0*px(3)/dd 662 px(1)=pw(2) 663 px(2)=pw(3) 664 px(3)=pw(1) 665 pn(1)=pw(2)*px(3)-pw(3)*px(2) 666 pn(2)=pw(3)*px(1)-pw(1)*px(3) 667 pn(3)=pw(1)*px(2)-pw(2)*px(1) 668 pz(1)=0.0d0 669 pz(2)=0.0d0 670 pz(3)=0.0d0 671 angle=0.95120d0 672 call rotate(pz,pn,angle,pw,dx) 673 xw(1,2,iwm)=xw(1,1,iwm)+dx(1) 674 xw(2,2,iwm)=xw(2,1,iwm)+dx(2) 675 xw(3,2,iwm)=xw(3,1,iwm)+dx(3) 676 angle=-0.95120d0 677 call rotate(pz,pn,angle,pw,dx) 678 xw(1,3,iwm)=xw(1,1,iwm)+dx(1) 679 xw(2,3,iwm)=xw(2,1,iwm)+dx(2) 680 xw(3,3,iwm)=xw(3,1,iwm)+dx(3) 681c 682 iwfnd(iwm)=3 683 endif 684 23 continue 685 endif 686c 687 if(util_print('restart',print_default)) then 688 write(lfnout,2003) 689 2003 format(' Crystal waters done',/) 690 endif 691c 692 argos_prep_misfit=.true. 693 return 694 9999 continue 695 argos_prep_misfit=.false. 696 return 697 end 698