PROGRAM DODCG1 IMPLICIT none integer*4 MX,iper,nt0,it,i,N,J,ii,nc,ns,jx,MXB,n6,npar0, 1ig,is,ix,NG,nat,ncs,ierr,ITER,NATB,nwl,nwl0,n4,n40, 2ivel,imag,ntr,nd PARAMETER (MX=4600,nt0=550,MXB=50,npar0=14,n40=6,nwl0=nt0) c Limits: c MX -# of atoms in peptide c nt0 -# of transitions on one chromophore c MXB -# of overlaped atoms c npar0 -# of parameters to follow c n40 -# of potential sites c nwl0 -# of fitted transitions integer*4 ic(MX),ics(MX),ttype(MX),IP(npar0) REAL*8 R(3*MX),rc(3),wc(MX),rn(3*MX),mag(nt0,3),elel(nt0,3), 1elev(nt0,3),en(nt0),en1(nt0,3),sum,w,A,XS,XB,RTOL,um(3),ue(3),dco, 1PG(nt0,npar0),zim(n40),dco0,eau, 1buvcd(nwl0,n40),fi(n40),F(npar0),enc(nt0), 1dcn,dcn0,dhn,dhn0,DPCX,DPCY,DPCZ,elelc(nt0,3),rcs(3) character*80 s80 character*1 prog logical lper,lfdp(4),lvel,llen,lene,lmag,lgeo COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NATB,ITER c c takes electronic dipoles from G.OUT, spreads them on large peptide c c Input: G.OUT -- TDDFT G Output c or moments -- TDDFT Turbomole output after transformation by program turbomoments c polymer.x --- the whole system c c Read in transition energies and dipoles from Gaussian or Turbomole: lvel=.false. llen=.false. lene=.false. lmag=.false. lgeo=.false. open(4,file='DODC.CON') read(4,400)prog 400 format(a1) if(prog.eq.'g'.or.prog.eq.'G')then open(30,file='G.OUT',status='old') 10 read(30,300,end=949,err=949)s80 300 format(a80) if(s80(12:44).eq.'excited state Transition electric')then call rd(nt0,it,elel,llen) write(6,*)it,' electric dipoles' goto 10 endif if(s80(12:44).eq.'excited state transition electric')then call rd(nt0,it,elel,llen) write(6,*)it,' electric dipoles' goto 10 endif if(s80(12:46).eq.'excited state transition velocity d')then call rd(nt0,it,elev,lvel) write(6,*)it,' electric velocity dipoles' goto 10 endif if(s80(12:44).eq.'excited state transition magnetic')then call rd(nt0,it,mag,lmag) write(6,*)it,' magnetic dipoles' goto 10 endif if(s80(2:15).eq.'Excited State '.and.s80(16:16).ne.'s')then lene=.true. read(s80(15:18),*)nd do 203 i=1,79 203 if(s80(i:i+1).eq.'nm')read(s80(i- 9:i-2),*)en(nd) endif if(s80(26:46).eq.'Standard orientation:')call ro(rn,nat,lgeo) goto 10 949 close(30) write(6,*)nd,' energy transitions in G.OUT' do 481 i=1,it c nm to hartrees: eau=1.0d7/(8065.54476345045d0*27.211384205943d0)/en(i) do 481 ix=1,3 c correct also gaussian sign: 481 elev(i,ix)=-elev(i,ix)/eau else if(prog.eq.'t'.or.prog.eq.'T')then open(30,file='moments',status='old') 100 read(30,300,end=9490,err=9490)s80 if(s80(2:46).eq.'Transition dipole moments in length formalism') 1 call rdt(nt0,it,elel,llen) if(s80(2:46).eq.'Transition dipole moments in velocity formali') 1 call rdt(nt0,it,elev,lvel) if(s80(2:28).eq.'Transition magnetic moments') 1 call rdt(nt0,it,mag,lmag) if(s80(2:24).eq.'Energies of Transitions') 1 call rdt(nt0,it,en1,lene) if(s80(2:19).eq.'Atomic coordinates')call rot(rn,nat,lgeo) goto 100 9490 close(30) do 9432 i=1,it 9432 en(i)=en1(i,2) write(6,*)it,' transitions in moments' else call report('unknown instruction '//prog) endif endif if(.not.llen)call report(' Electric length dipoles not found') if(.not.lvel)call report(' Electric velocity dipoles not found') if(.not.lmag)call report(' Magnetic dipoles not found') if(.not.lene)call report(' Energies not found') if(.not.lgeo)call report(' Geometry not found') c read-in geometry OPEN(2,FILE='BIG.X',STATUS='OLD') read(2,*) read(2,*)N if(N.gt.MX)call report('Too many atoms.') do 1 i=1,N 1 READ(2,*)ttype(i),(R(J+3*(i-1)),J=1,3) close(2) write(6,*)N, ' atoms in BIG.X' c read(4,*)ivel,imag if(ivel.eq.1)then write(6,*)' velocity dipoles used' do 48 i=1,it do 48 ix=1,3 48 elel(i,ix)=elev(i,ix) else write(6,*)' length dipoles used' endif c c transform aus to Bohr magnetons, debyes do 49 i=1,it do 49 ix=1,3 c this is not sure, check the magnetic moment units!: mag(i,ix)=mag(i,ix)/2.0d0*2.541765d0 49 elel(i,ix)=elel(i,ix)*2.541765d0 if(imag.eq.1)then write(6,*)' magnetic moments added' else write(6,*)' magnetic moments skipped' do 47 i=1,it do 47 ix=1,3 47 mag(i,ix)=0.0d0 endif open(3,file='DC.SC') c chromophore atoms: H N C O read(4,*)ncs,(ics(i),i=1,ncs) dco0=dsqrt((rn(1+3*(ics(4)-1))-rn(1+3*(ics(3)-1)))**2+ 1 (rn(2+3*(ics(4)-1))-rn(2+3*(ics(3)-1)))**2+ 1 (rn(3+3*(ics(4)-1))-rn(3+3*(ics(3)-1)))**2) dcn0=dsqrt((rn(1+3*(ics(2)-1))-rn(1+3*(ics(3)-1)))**2+ 1 (rn(2+3*(ics(2)-1))-rn(2+3*(ics(3)-1)))**2+ 1 (rn(3+3*(ics(2)-1))-rn(3+3*(ics(3)-1)))**2) dhn0=dsqrt((rn(1+3*(ics(2)-1))-rn(1+3*(ics(1)-1)))**2+ 1 (rn(2+3*(ics(2)-1))-rn(2+3*(ics(1)-1)))**2+ 1 (rn(3+3*(ics(2)-1))-rn(3+3*(ics(1)-1)))**2) read(4,*)NG write(3,*)NG,' groups' write(6,*)NG,' groups' c c geometry and ES field corrections: read(4,*)(IP(i),i=1,npar0) write(6,6009)(IP(i),i=1,npar0) 6009 format(14i2) read(4,*)ntr do 6 it=1,ntr 6 read(4,*)(PG(it,i),i=1,npar0) c c determine which reference value to use if(IP(1).lt.0)then write(6,*)' Ab initio value of d(C=O):' else dco0=1.231d0 write(6,*)' MD value of d(C=O):' endif write(6,6070)dco0 6070 format(f12.3,' A') if(IP(3).lt.0)then write(6,*)' Ab initio value of d(C=N):' else dcn0=1.335d0 write(6,*)' MD value of d(C=N):' endif write(6,6070)dcn0 c c periodic boundary conditions read(4,*)iper,DPCX,DPCY,DPCZ lper=iper.eq.1 c c potential correction: if(IP(npar0).eq.1)then open(55,file='BIJ.SCR',status='old',form='unformatted') read(55)n4,n6,(zim(i),i=1,n4),(lfdp(i),i=1,4) write(6,*)'Zs:' write(6,6071)(zim(i),i=1,n4) 6071 format(6f12.4) if(lfdp(1))read(55) if(lfdp(2))read(55) if(lfdp(3))read(55) if(lfdp(4))read(55)nwl,((buvcd(i,j),i=1,nwl),j=1,n6) close(55) endif do 2 IG=1,NG write(3,*)IG,' . group:' write(6,*)IG,' . group:' c Number of atoms determining center and their numbers: read(4,*)nc,(ic(i),i=1,nc) write(6,*)nc,(ic(i),i=1,nc) if(nc.ne.ncs)call report('nc<>ncs!') c center weights: read(4,*)(wc(i),i=1,nc) sum=0.0d0 do 321 ii=1,nc 321 sum=sum+wc(ii) do 3 i=1,3 rc(i)=0.0d0 rcs(i)=0.0d0 do 32 ii=1,nc rcs(i)=rcs(i)+rn(i+3*(ics(ii)-1))*wc(ii) 32 rc(i)=rc(i)+R(i+3*(ic(ii)-1))*wc(ii) rcs(i)=rcs(i)/sum 3 rc(i)=rc(i)/sum write(3,1000)(rc(i)/10.0d0,i=1,3) 1000 format(3f15.6,' nm') ns=ntr+1 write(3,*)ns,' states' write(3,*) c c ground state: write(3,*)'0 0.0' c c correct energies and dipoles according to c geometry and ES field:w=w0*(1+dd/d*p) dip=dip0*(1+dd/d*p) c first par-energy, second dipole c do 41 is=2,ns do 46 i=1,npar0 46 F(i)=0.0d0 if(IP(1).ne.0)then dco=dsqrt((R(1+3*(ic(4)-1))-R(1+3*(ic(3)-1)))**2+ 1 (R(2+3*(ic(4)-1))-R(2+3*(ic(3)-1)))**2+ 1 (R(3+3*(ic(4)-1))-R(3+3*(ic(3)-1)))**2) F(1)=(dco-dco0)/dco0*PG(is-1,1) F(2)=(dco-dco0)/dco0*PG(is-1,2) if(is.eq.2)write(6,6001)1,dco,dco0 6001 format(i6,' d(C=O) = ',f10.5,'; D0 = ',f10.5) endif c C=N bond if(IP(3).ne.0)then dcn=dsqrt((R(1+3*(ic(2)-1))-R(1+3*(ic(3)-1)))**2+ 1 (R(2+3*(ic(2)-1))-R(2+3*(ic(3)-1)))**2+ 1 (R(3+3*(ic(2)-1))-R(3+3*(ic(3)-1)))**2) F(3)=(dcn-dcn0)/dcn0*PG(is-1,3) F(4)=(dcn-dcn0)/dcn0*PG(is-1,4) if(is.eq.2)write(6,6003)3,dcn,dcn0 6003 format(i6,' d(C-N) = ',f10.5,'; D0 = ',f10.5) endif c H-N bond if(IP(5).ne.0)then dhn=dsqrt((R(1+3*(ic(2)-1))-R(1+3*(ic(1)-1)))**2+ 1 (R(2+3*(ic(2)-1))-R(2+3*(ic(1)-1)))**2+ 1 (R(3+3*(ic(2)-1))-R(3+3*(ic(1)-1)))**2) F(5)=(dhn-dhn0)/dhn0*PG(is-1,5) F(6)=(dhn-dhn0)/dhn0*PG(is-1,6) if(is.eq.2)write(6,6004)5,dhn,dhn0 6004 format(i6,' d(N-H) = ',f10.5,'; D0 = ',f10.5) endif enc(is-1)=en(is-1)*(1.0d0+F(1)+F(3)+F(5)) do 41 ix=1,3 41 elelc(is-1,ix)=elel(is-1,ix)*(1.0d0+F(2)+F(4)+F(6)) if(IP(npar0).ne.1)then call fw(n6,fi,N,ttype,R,ic(1),ic(2),ic(3),ic(4),lper, 1 DPCX,DPCY,DPCZ) do 42 is=2,ns if(en(is-1).gt.200.0d0)then do 43 i=1,4 43 enc(is-1)=enc(is-1)+buvcd(1,i)*fi(i) else do 44 i=1,4 44 enc(is-1)=enc(is-1)+buvcd(2,i)*fi(i) endif 42 continue endif do 4 is=2,ns w=1.0d7/enc(is-1) 4 write(3,1001)is,w/8065.54093D0,w,enc(is-1) 1001 format(i3,f11.6,' eV =',f11.2,' cm-1 =',f11.2,' nm') write(3,*) c c write group reference atoms into XS and corresponding ones in XB: if(ncs.gt.MXB)call report('ncs>MXB') do 5 i=1,ncs do 51 ix=1,3 XB(i,ix)= R(ix+3*(ic(i)-1))-rc(ix) 51 XS(i,ix)=rn(ix+3*(ics(i)-1))-rcs(ix) 5 continue NATB=ncs c c overlap XS,XB, make the rotation matrix, store in A: call domatrix(ierr) if(ierr.ne.0)call report('rotation not successful') write(6,*)' Transformation matrix' do 67671 i=1,3 67671 write(6,6767)(A(i,j),j=1,3) 6767 format(3f10.5) c c write rotated dipole moments: do 2 is=2,ns do 8 ix=1,3 ue(ix)=0.0d0 um(ix)=0.0d0 do 8 jx=1,3 ue(ix)=ue(ix)+A(ix,jx)*elelc(is-1,jx) 8 um(ix)=um(ix)+A(ix,jx)* mag(is-1,jx) c 2 write(3,1002)is,(ue(ix),ix=1,3),(um(ix),ix=1,3) 1002 format(i3,6f12.6) close(3) close(4) call report('DC.SC written') end c SUBROUTINE REPORT(RS) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*(*) RS WRITE(6,3000) 3000 FORMAT(/,80(1H*)) WRITE(6,*)RS WRITE(6,3001) 3001 FORMAT(80(1H*),/,/,'PROGRAM STOPPED') STOP END subroutine rd(nt0,it,e,l) implicit none integer*4 nt0,it,idumm,i character*80 s80 real*8 e(nt0,3) logical nu,l l=.true. it=0 read(30,*) 31 read(30,300)s80 300 format(a80) if(nu(s80(10:10)))then it=it+1 if(it.gt.nt0)call report(' too many transitions !') read(s80,*)idumm,(e(it,i),i=1,3) goto 31 else backspace 30 return endif end subroutine rdt(nt0,it,e,l) implicit none integer*4 nt0,it,idumm,i character*10 s10 real*8 e(nt0,3) logical nu,l l=.true. it=0 read(30,*) 31 read(30,300)s10 300 format(a10) if(nu(s10(4:4)))then it=it+1 if(it.gt.nt0)then write(6,*)' too many transitions !' stop endif backspace 30 read(30,*)idumm,(e(it,i),i=1,3) goto 31 endif backspace 30 return end subroutine ro(rn,nat,l) implicit none integer*4 dum,i,nat character*80 s80 real*8 rn(*) logical nu,l l=.true. nat=0 read(30,*) read(30,*) read(30,*) read(30,*) 31 read(30,300)s80 300 format(a80) if(nu(s80(5:5)).or.nu(s80(7:7)))then nat=nat+1 read(s80,*)dum,dum,dum,(rn(i+3*(nat-1)),i=1,3) goto 31 endif write(6,*)nat,' atoms in chromophore' return end subroutine rot(rn,nat,l) implicit none integer*4 i,nat character*80 s80 character*15 dum real*8 rn(*) logical l l=.true. nat=0 read(30,*) 31 read(30,300)s80 300 format(a80) if(s80(1:15).ne.' ')then nat=nat+1 read(s80,*)dum,(rn(i+3*(nat-1)),i=1,3) goto 31 endif write(6,*)nat,' atoms in chromophore' return end subroutine re(e) implicit none integer*4 it character*80 s80 real*8 e(*) backspace 30 read(30,300)s80 300 format(a80) read(s80(15:18),*)it read(s80(47:55),*)e(it) return end function nu(c) character*1 c logical nu nu=c.eq.'1'.or.c.eq.'2'.or.c.eq.'3'.or.c.eq.'4'.or.c.eq.'5'. 1or.c.eq.'6'.or.c.eq.'7'.or.c.eq.'8'.or.c.eq.'9'.or.c.eq.'0' return end SUBROUTINE DOMATRIX(IERR) C This is a Fortran version of Petr Malon's subroutine Amoeba IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MXB=50) DIMENSION ANG(3),P(4,3),Y(4),PBAR(3),PR(3),PRR(3) COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER IERR=0 C STARTING VALUES FOR THE ITERATION: ANG(1)=0.2d0 ANG(2)=0.0d0 ANG(3)=0.0d0 DO 1 I=1,3 1 P(1,I)=ANG(I) Y(1)=FU(ANG) ANG(2)=ANG(1) ANG(1)=0.0d0 DO 2 I=1,3 2 P(2,I)=ANG(I) Y(2)=FU(ANG) ANG(3)=ANG(2) ANG(2)=0.0d0 DO 3 I=1,3 3 P(3,I)=ANG(I) Y(3)=FU(ANG) ANG(1)=ANG(3) ANG(2)=ANG(3) ANG(3)=0.0d0 DO 4 I=1,3 4 P(4,I)=ANG(I) Y(4)=FU(ANG) FTOL=0.0000001d0 ITMAX=1500 ITER=0 99999 ILO=1 IF(Y(1).GT.Y(2))THEN IHI=1 INHI=2 ELSE IHI=2 INHI=1 ENDIF DO 5 I=1,4 IF(Y(I).LT.Y(ILO))ILO=I IF(Y(I).GT.Y(IHI))THEN INHI=IHI IHI=I ELSE IF(Y(I).GT.Y(INHI))THEN IF(I.NE.IHI)INHI=I ENDIF ENDIF 5 CONTINUE IF(ABS(Y(IHI))+ABS(Y(ILO)).LT.0.00000000001d0)THEN RTOL=0.0d0 ELSE RTOL=ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))*2.0d0 ENDIF IF(RTOL.LT.FTOL)RETURN IF(ITER.EQ.ITMAX)THEN WRITE(3,*)' Rotation has not converged !' IERR=1 RETURN ENDIF ITER=ITER+1 DO 55 I=1,3 55 PBAR(I)=0.0d0 DO 6 I=1,4 IF(I.NE.IHI)THEN DO 7 J=1,3 7 PBAR(J)=PBAR(J)+P(I,J) ENDIF 6 CONTINUE DO 8 J=1,3 PBAR(J)=PBAR(J)/3.0d0 8 PR(J)=2.0d0*PBAR(J)-P(IHI,J) YPR=FU(PR) IF(YPR.LE.Y(ILO))THEN DO 9 J=1,3 9 PRR(J)=2.0d0*PR(J)-PBAR(J) YPRR=FU(PRR) IF(YPRR.LT.Y(ILO))THEN DO 10 J=1,3 10 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 11 J=1,3 11 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ELSE IF(YPR.GE.Y(INHI))THEN IF(YPR.LT.Y(IHI))THEN DO 12 J=1,3 12 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF DO 13 J=1,3 13 PRR(J)=0.5d0*P(IHI,J)+0.5d0*PBAR(J) YPRR=FU(PRR) IF(YPRR.LT.Y(IHI))THEN DO 14 J=1,3 14 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 15 I=1,4 IF(I.NE.ILO)THEN DO 16 J=1,3 PR(J)=0.5d0*(P(I,J)+P(ILO,J)) 16 P(I,J)=PR(J) Y(I)=FU(PR) ENDIF 15 CONTINUE ENDIF ELSE DO 17 J=1,3 17 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ENDIF GOTO 99999 END FUNCTION FU(ANG) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MXB=50) COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER DIMENSION ANG(3) S1=SIN(ANG(1)) S2=SIN(ANG(2)) S3=SIN(ANG(3)) C1=COS(ANG(1)) C2=COS(ANG(2)) C3=COS(ANG(3)) C ang1 = theta C ang2 = phi C ang3 = cappa A(1,1)=C1*C2*C3-S2*S3 A(1,2)=C1*S2*C3+C2*S3 A(1,3)=-S1*C3 A(2,1)=-C1*C2*S3-S2*C3 A(2,2)=-C1*S2*S3+C2*C3 A(2,3)=S1*S3 A(3,1)=S1*C2 A(3,2)=S1*S2 A(3,3)=C1 R=0.0d0 DO 1 I=1,NAT TX=A(1,1)*XS(I,1)+A(1,2)*XS(I,2)+A(1,3)*XS(I,3) TY=A(2,1)*XS(I,1)+A(2,2)*XS(I,2)+A(2,3)*XS(I,3) TZ=A(3,1)*XS(I,1)+A(3,2)*XS(I,2)+A(3,3)*XS(I,3) 1 R=R+(TX-XB(I,1))*(TX-XB(I,1))+ 1 (TY-XB(I,2))*(TY-XB(I,2))+(TZ-XB(I,3))*(TZ-XB(I,3)) FU=R RETURN END c subroutine fw(n6,fi,nat,ity,r,iih,iin,iic,iio,lper, 1DPCX,DPCY,DPCZ) c find waters and their potentials implicit none integer*4 i,n6,jj,ity(*),ityo,ityh,j1,nh,kk,k1,iih,iin,iic,iio, 1ii,nat logical lper real*8 fi(*),xo,yo,zo,r(*),xh,yh,zh,d2,x1,y1,z1, 1x2,y2,z2,qh,qo,DPCX,DPCY,DPCZ,dx2,dy2,dz2 data ityo,ityh,qh,qo/649,650,0.41d0,-0.82d0/ dx2=DPCX/2.0d0 dy2=DPCY/2.0d0 dz2=DPCZ/2.0d0 do 4 i=1,n6 4 fi(i)=0.0d0 do 6 jj=1,nat if(ity(jj).eq.ityo)then j1=3*(jj-1) xo=r(j1+1) yo=r(j1+2) zo=r(j1+3) nh=0 do 7 kk=1,nat if(ity(kk).eq.ityh)then k1=3*(kk-1) xh=r(k1+1) yh=r(k1+2) zh=r(k1+3) d2=(xo-xh)**2+(yo-yh)**2+(zo-zh)**2 if(d2.gt.0.64d0.and.d2.lt.1.44d0)then nh=nh+1 if(nh.gt.2)call report('Too many Hs around one O') if(nh.eq.1)then x1=xh y1=yh z1=zh else x2=xh y2=yh z2=zh endif endif endif 7 continue if(nh.eq.2)then call add(iih ,1,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) call add(iin ,2,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) call add(iic ,3,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) call add(iio ,4,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) c call add(iian,5,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, c 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) c call add(iiac,6,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, c 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) endif endif 6 continue write(6,*)' Potentials:' write(6,1603)(fi(ii),ii=1,n6) 1603 format(8f9.3) return end c subroutine add(ii,ic,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) implicit none real*8 r(*),xo,yo,zo,x1,y1,z1,x2,y2,z2,fi(*),x,y,z,ro,r1,r2,qh,qo, 1dx2,dy2,dz2,DPCX,DPCY,DPCZ integer*4 ii,ic logical lper c x=r(3*(ii-1)+1) y=r(3*(ii-1)+2) z=r(3*(ii-1)+3) if(lper)then if(x-xo.gt.dx2.and.x-x1.gt.dx2.and.x-x2.gt.dx2)x=x-DPCX if(y-yo.gt.dy2.and.y-y1.gt.dy2.and.y-y2.gt.dy2)y=y-DPCY if(z-zo.gt.dz2.and.z-z1.gt.dz2.and.z-z2.gt.dz2)z=z-DPCZ if(xo-x.gt.dx2.and.x1-x.gt.dx2.and.x2-x.gt.dx2)x=x+DPCX if(yo-y.gt.dy2.and.y1-y.gt.dy2.and.y2-y.gt.dy2)y=y+DPCY if(zo-z.gt.dz2.and.z1-z.gt.dz2.and.z2-z.gt.dz2)z=z+DPCZ endif ro=dsqrt((xo-x)**2+(yo-y)**2+(zo-z)**2) r1=dsqrt((x1-x)**2+(y1-y)**2+(z1-z)**2) r2=dsqrt((x2-x)**2+(y2-y)**2+(z2-z)**2) c fi(ic)=fi(ic)+qh/r1+qh/r2+qo/ro return end