C C $Id: atmchg.F,v 1.5 1998/08/14 17:49:01 arjan Exp arjan $ C C------------------------------------------------------------------------ subroutine pme_calccm1 implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" #ifdef MPI #include "mpif.h" #endif dimension sp2(maxatm),bksave(maxatm),dqsave(maxatm),cksave(maxatm) logical calccm1, calccm2 CC-----------------------------------------------------------CC IIMAX = IIMAT(NATOMS+1)-1 IPMAX = IP1(NATOMS+1) IJMAX = IJMAT(IPMAX)-1 do i=1,mxpair bkk1(i)=0.0d0 enddo DO II=1,IIMAX FPMEDIAG(II) = 0.0D0 ENDDO DO IJ=1,IJMAX FPMEDIAT1(IJ) = 0.0D0 FPMEDIAT2(IJ) = 0.0D0 ENDDO calccm1 = .true. calccm2 = .false. C get the Mulliken charges ii = 0 do 20 i=1,natoms iai = iatnum(i) c skip dummy atom. if(iai.eq.0) go to 20 norbsi = natorb(iai) zi = zchg(iai) psumii = 0.0d0 if(norbsi.gt.0)then do 10 iiorb=1,norbsi ii = ii + iiorb psumii = psumii + pdiag(ii) 10 continue endif atchg(i) = zi - psumii atchg2(i) = atchg(i) 20 continue C for now atchg refers to Mulliken charges, atchg2 to CM1 charges, C atchg3 to CM2 charges (swap this if needed at the end of the C routine) C . there are two efficient ways to calculate the CM1 charges: C . the faster one is a loop from iat=2,natoms with C . an inner loop from jat=1,iat-1, calculate and store the C . B(iat,jat) = B(jat, iat) coefficients and then C . do another iat=2,natoms - jat=1,iat-1 loop for C . the normalization using the stored B(iat,jat) = B(jat,iat) C . coefficients. C . The second slower one loops from iat=1,natoms in which the C . inner loop jat=1,natoms calculates the B(iat,jat) coefficients. C . After this jat-loop a partial normalization is performed. C . Method one is faster because it only performs the C . density matrix summation needed for the B(iat,jat) calculation C . N^2/2 times (method 2 does this N^2 times), C . but requires a lot more memory: N^2/2 elements to C . store the B(iat,jat) coefficients. Method 2 only needs to C . store the N intermediate B(iat,jat) coefficients within the C . iat-loop (N is the number of atoms). C . Here memory considerations lead me to implement the second C . (slower) method. C . CM2 charges can be cheaply obtained in the same sweep. C . calculate CM1P / CM2P charges do 2000 iat=1,natoms iatk = iatnum(iat) if (iatk.eq.0) goto 2000 bk = 0.0 bhk = 0.0 norbsi = natorb(iatk) npairs = ip1(natoms+1)-1 do 200 jat=1,natoms sp2(jat) = 0.0 if (jat.eq.iat) goto 200 jatk = iatnum(jat) if (jatk.ne.0) then c i = max(iat,jat) c j = min(iat,jat) if (iat.gt.jat) then call ijfind(npairs,iat,jat,ijaddr) else call ijfind(npairs,jat,iat,ijaddr) endif if (ijaddr.ne.0) then norbsj = natorb(jatk) norbsij = norbsi*norbsj - 1 ijstart = ijmat(ijaddr) ijend = ijstart + norbsij do ij=ijstart,ijend p2 = pdiat(ij)*pdiat(ij) sp2(jat) = sp2(jat) + p2 bk = bk + p2 enddo call pme_getcm1p(iatk,jat,jatk,bhk,sp2) endif endif 200 continue call pme_calccm1p(iat,iatk,bhk,bk,sp2,ck,dq) dqsave(iat)=dq bksave(iat)=bk cksave(iat)=ck c write(59,*)dqsave(iat),bksave(iat),cksave(iat) 2000 continue #ifdef MPI DO IAT=myid+1,NATOMS,nproc #else DO IAT=1,NATOMS #endif ! do iat=1,natoms iatk = iatnum(iat) if (iatk.ne.0) then norbsi = natorb(iatk) IJDIAG = IIMAT(IAT)-1 c write(59,*)iat,jat,bk do iorb=1,norbsi ijdiag=ijdiag+iorb FPMEDIAG(IJDIAG)= -1 -bksave(iat)*cksave(iat) enddo endif enddo ijbond=1 do 2001 iat=2,natoms iatk = iatnum(iat) if (iatk.eq.0) goto 2001 norbsi = natorb(iatk) do 201 jat=1,iat-1 c write(*,*)iat,jat c if (jat.eq.iat) goto 201 jatk = iatnum(jat) if (jatk.ne.0) then c i = max(iat,jat) c j = min(iat,jat) call ijfind(npairs,iat,jat,ijaddr) if (ijaddr.ne.0) then c call ijfind(npairs,j,i,jiaddr) norbsj = natorb(jatk) c norbsij = norbsi*norbsj - 1 c ijstart = ijmat(ijaddr) c ijend = ijstart + norbsij c jistart = ijmat(jiaddr) c jiend = jistart + norbsij IJ1= IJMAT(ijaddr) IJ=IJMAT(ijaddr) call pme_getdhk(iatk,jatk,dhkj) call pme_getdhk(jatk,iatk,dhki) bkk1(ijaddr)=0.0d0 do iorbi=1,norbsi do iorbj=1,norbsj c do iorb=ijstart,ijend bkk1(ijaddr) = bkk1(ijaddr) + pdiat(ij1)*pdiat(ij1) c enddo ij1=ij1+1 enddo enddo do iorbi=1,norbsi do iorbj=1,norbsj fpmediat1(ij) = bksave(iat)*dhkj*pdiat(ij) & - bkk1(ijaddr)*dhki*pdiat(ij) & - pdiat(ij)*dqsave(jat) & + pdiat(ij)*dqsave(iat) c fpmediat2(ij) = -dhkj*bkk*2.0d0*pdiat(ij) c & -2.0d0*pdiat(ij)*dqsave(iat)-2.0d0*pdiat(ij)*dqsave(iat) fpmediat2(ij) = bksave(jat)*dhki*pdiat(ij) & - bkk1(ijaddr)*dhkj*pdiat(ij) & - pdiat(ij)*dqsave(iat) & + pdiat(ij)*dqsave(jat) c fpmediat2(ij) = -fpmediat1(ij) c fpmediat2(ij) = dhki*bkk*pdiat(ij) c & - pdiat(ij)*dqsave(iat)+ pdiat(ij)*dqsave(jat) c & + 2.0d0*pdiat(ij)*dqsave(i) ij=ij+1 c write(59,*) fpmediat1(ij),pdiat(ij),dqsave(i) enddo enddo c call pme_getdhk(iatk,jatk,dhk) c do ij=ijstart,ijend c fpmediat(ij) = bksave(iat)*dhk*2.0d0*pdiat(ij) c & + dqsave(iat)*2.0d0*pdiat(ij) - 2.0d0*pdiat(ij)*dqsave(jat) c write(59,*) ij, fpmediat(ij) c enddo c do ji=jistart,jiend c fpmediat(ji) = fpmediat(ji) - 2.0d0*pdiat(ji)*dq c enddo endif endif ijbond=ijbond+1 201 continue 2001 continue c write(59,*) (fpmediat2(i),i=1,100) end CC--------------------------------------------------------------------CC CC--------------------------------------------------------------------CC C C subroutines for CM1/PM3 charge model C CC--------------------------------------------------------------------CC CC--------------------------------------------------------------------CC subroutine pme_getcm1p(iatk,jat,jatk,bhk,sp2) implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" dimension sp2(maxatm),dhk(maxatm) if (iatk.eq.1) then C . H atom if (jatk.eq.7) then C . overlap with N atom d = 0.1854d0 elseif (jatk.eq.8) then C . overlap with O atom d = 0.1434d0 elseif (jatk.eq.14) then C . overlap with Si atom d = -0.1004d0 else d = 0.0d0 endif c d=0.0d0 bhk = bhk + d*sp2(jat) endif end CC--------------------------------------------------------------------CC subroutine pme_calccm1p(iat,iatk,bhk,bk,sp2,c,dq) implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" dimension sp2(maxatm) if (iatk.eq.7) then C . N atom c = 0.0d0 d = -0.0909d0 elseif (iatk.eq.8) then C . O atom c = 0.0d0 d = -0.0449d0 elseif (iatk.eq.9) then C . F atom c = 0.3381d0 d = 0.0148d0 elseif (iatk.eq.16) then C . S atom c = -0.0834d0 d = -0.0848d0 elseif (iatk.eq.17) then C . Cl atom c = -0.1080d0 d = -0.1168d0 elseif (iatk.eq.35) then C . Br atom c = -0.0116d0 d = -0.0338d0 elseif (iatk.eq.53) then C . I atom c = -0.3213d0 d = -0.0636d0 else c = 0.0d0 d = 0.0d0 endif c d=0.0d0 c c=0.0d0 dq = c*atchg(iat) + d + bhk atchg2(iat) = atchg2(iat) + dq*bk C partial normalization: do jat=1,natoms atchg2(jat) = atchg2(jat) - sp2(jat)*dq enddo end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine pme_getdhk(iatk,jatk,dhk) implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" if (iatk.eq.1) then C . H atom if (jatk.eq.7) then C . overlap with N atom d = 0.1854d0 elseif (jatk.eq.8) then C . overlap with O atom d = 0.1434d0 elseif (jatk.eq.14) then C . overlap with Si atom d = -0.1004d0 else d = 0.0d0 endif else d=0.0d0 endif c d=0.0d0 dhk=d end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC