C C $Id: atmchg.F,v 1.5 1998/08/14 17:49:01 arjan Exp arjan $ C C------------------------------------------------------------------------ subroutine atmchg implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" dimension sp2(maxatm) logical calccm1, calccm2 CC-----------------------------------------------------------CC C for non-MC simulations, CM1, CM2 and Mulliken charges C are calculated. C atchg is the charge array used in the PME calculations. C so it depends on the value of cm1 and cm2 if atchg will be filled C with Mulliken charges or CM1 or CM2 charges. C The following table shows which charges will be assigned to C which array (* is wildcard, valid for both T(rue) and F(false), C - is not calculated) C C no.|mcsim | pme | cm1 | cm2 | Mulliken | CM1 | CM2 C 1 T * F F atchg - - C 2 T * T F - atchg - C 3 T * F T - - atchg C C 4 F F * * atchg atchg2 atchg3 C 5 F T F F atchg atchg2 atchg3 C 6 F T T F atchg2 atchg atchg3 C 7 F T F T atchg3 atchg2 atchg C if (.not.mcsim) then calccm1 = .true. calccm2 = .true. elseif (cm1) then calccm1 = .true. calccm2 = .false. else if (cm2) then calccm1 = .false. calccm2 = .true. else calccm1 = .false. calccm2 = .false. endif 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) atchg3(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) if (calccm1.or.calccm2) then 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. if (am1) then C . calculate the CM1A / CM2A charges do 1000 iat=1,natoms iatk = iatnum(iat) if (iatk.eq.0) goto 1000 bk = 0.0 bnc1 = 0.0 bnc2 = 0.0 bos = 0.0 bno = 0.0 bhk = 0.0 norbsi = natorb(iatk) npairs = ip1(natoms+1)-1 do 100 jat=1,natoms sp2(jat) = 0.0 if (jat.eq.iat) goto 100 jatk = iatnum(jat) if (jatk.ne.0) then i = max(iat,jat) j = min(iat,jat) call ijfind(npairs,i,j,ijaddr) 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 if (calccm1) & call getcm1a(iatk,jat,jatk,bhk, & bnc1,bnc2,bno,bos,sp2) if (calccm2) & call getcm2a(iat,iatk,jat,jatk,sp2) endif endif 100 continue if (calccm1) & call calccm1a(iat,iatk,bhk,bnc1,bnc2, & bno,bos,bk,sp2) 1000 continue elseif (pm3) then 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 i = max(iat,jat) j = min(iat,jat) call ijfind(npairs,i,j,ijaddr) 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 if (calccm1) & call getcm1p(iatk,jat,jatk,bhk,sp2) if (calccm2) & call getcm2p(iat,iatk,jat,jatk,sp2) endif endif 200 continue if (calccm1) & call calccm1p(iat,iatk,bhk,bk,sp2) 2000 continue endif C . see if swapping of the charge arrays is needed C . right now atchg contains the Mulliken charges, atchg2 the CM1 C . and atchg3 the CM2 charges C . see table in beginning of subroutine C . no. 1, 4 and 5 don't require any swapping if (cm1) then C . no. 2 and 6 C . swap atchg with atchg2 do i=1,natoms tmp = atchg(i) atchg(i) = atchg2(i) atchg2(i) = tmp enddo elseif (cm2) then C . no. 3 and 7 C . swap atchg with atchg3 do i=1,natoms tmp = atchg(i) atchg(i) = atchg3(i) atchg3(i) = tmp enddo endif endif end CC--------------------------------------------------------------------CC CC--------------------------------------------------------------------CC C C subroutines for CM1/AM1 charge model C CC--------------------------------------------------------------------CC CC--------------------------------------------------------------------CC subroutine getcm1a(iatk,jat,jatk,bhk,bnc1,bnc2,bno,bos,sp2) implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" dimension sp2(maxatm) if (iatk.eq.1) then C . H atom if (jatk.eq.7) then C . overlap with N atom d = 0.0850 elseif (jatk.eq.8) then C . overlap with O atom d = 0.1447 elseif (jatk.eq.14) then C . overlap with Si atom d = 0.0640 else d = 0.0 endif bhk = bhk + d*sp2(jat) elseif (iatk.eq.7) then C . N atom if (jatk.eq.6) then C . overlap with C atom t = 10.0*(sp2(jat)-2.3) t = 0.5*tanh(t)+0.5 bnc1 = bnc1 + t bnc2 = bnc2 - 0.0880*t elseif (jatk.eq.8) then C . overlap with O atom bno = bno - 0.0630*sp2(jat) endif elseif ((iatk.eq.8).and.(jatk.eq.16)) then C . overlap S atom with O atom bos = bos - 0.060*sp2(jat) endif end CC--------------------------------------------------------------------CC subroutine calccm1a(iat,iatk,bhk,bnc1,bnc2,bno,bos,bk,sp2) implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" dimension sp2(maxatm) bnc1 = -0.3846*bnc1 if (iatk.eq.7) then C . N atom c = 0.3846 d = 0.0 elseif (iatk.eq.8) then C . O atom c = 0.0 d = -0.0283 elseif (iatk.eq.9) then C . F atom c = 0.1468 d = 0.0399 elseif (iatk.eq.16) then C . S atom c = -0.1311 d = -0.0956 elseif (iatk.eq.17) then C . Cl atom c = 0.0405 d = -0.0276 elseif (iatk.eq.35) then C . Br atom c = 0.1761 d = -0.0802 elseif (iatk.eq.53) then C . I atom c = 0.2380 d = -0.1819 else c = 0.0 d = 0.0 endif dq = (c + bnc1)*atchg(iat) + d + bhk + bos + & bnc2 + bno atchg2(iat) = atchg2(iat) + dq*bk C partial normalization: do jat=1,natoms atchg2(jat) = atchg2(jat) - sp2(jat)*dq enddo end CC--------------------------------------------------------------------CC CC--------------------------------------------------------------------CC C C subroutines for CM1/PM3 charge model C CC--------------------------------------------------------------------CC CC--------------------------------------------------------------------CC subroutine getcm1p(iatk,jat,jatk,bhk,sp2) implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" dimension sp2(maxatm) if (iatk.eq.1) then C . H atom if (jatk.eq.7) then C . overlap with N atom d = 0.1854 elseif (jatk.eq.8) then C . overlap with O atom d = 0.1434 elseif (jatk.eq.14) then C . overlap with Si atom d = -0.1004 else d = 0.0 endif bhk = bhk + d*sp2(jat) endif end CC--------------------------------------------------------------------CC subroutine calccm1p(iat,iatk,bhk,bk,sp2) 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.0 d = -0.0909 elseif (iatk.eq.8) then C . O atom c = 0.0 d = -0.0449 elseif (iatk.eq.9) then C . F atom c = 0.3381 d = 0.0148 elseif (iatk.eq.16) then C . S atom c = -0.0834 d = -0.0848 elseif (iatk.eq.17) then C . Cl atom c = -0.1080 d = -0.1168 elseif (iatk.eq.35) then C . Br atom c = -0.0116 d = -0.0338 elseif (iatk.eq.53) then C . I atom c = -0.3213 d = -0.0636 else c = 0.0 d = 0.0 endif 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 CC--------------------------------------------------------------------CC CC--------------------------------------------------------------------CC C C subroutines for CM2/AM1 charge model C CC--------------------------------------------------------------------CC CC--------------------------------------------------------------------CC subroutine getcm2a(iat,iatk,jat,jatk,sp2) implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" dimension sp2(maxatm) i = min(iatk,jatk) j = max(iatk,jatk) c = 0.0 d = 0.0 if (i.eq.1) then C . H atom if (j.eq.6) then C . overlap with C atom c = -0.020 elseif (j.eq.7) then C . overlap with N atom c = 0.207 elseif (j.eq.8) then C . overlap with O atom c = 0.177 elseif (j.eq.14) then C . overlap with Si atom c = -0.083 elseif (j.eq.15) then C . overlap with P atom d = 0.103 elseif (j.eq.16) then C . overlap with S atom c = 0.038 endif elseif (i.eq.6) then C . C atom if (j.eq.7) then C . overlap with N atom c = 0.008 d = 0.086 elseif (j.eq.8) then C . overlap with O atom c = 0.026 d = 0.016 elseif (j.eq.9) then C . overlap with F atom d = 0.019 elseif (j.eq.14) then C . overlap with Si atom c = 0.062 elseif (j.eq.15) then C . overlap with P atom d = -0.019 elseif (j.eq.16) then C . overlap with S atom c = -0.059 d = 0.171 elseif (j.eq.17) then C . overlap with Cl atom d = 0.027 elseif (j.eq.35) then C . overlap with Br atom d = 0.081 elseif (j.eq.53) then C . overlap with I atom d = 0.147 endif elseif (i.eq.7) then C . N atom if (j.eq.8) then C . overlap with O atom c = -0.197 d = 0.134 endif elseif (i.eq.8) then C . O atom if (j.eq.15) then C . overlap with P atom d = 0.088 endif elseif (i.eq.9) then C . F atom if (j.eq.15) then C . overlap with P atom d = 0.252 endif elseif (i.eq.15) then C . P atom if (j.eq.16) then C . overlap with S atom d = -0.080 endif endif t = sp2(jat)*(d + c*sp2(jat)) if (i.ne.iatk) t = -t atchg3(iat) = atchg3(iat) + t end CC--------------------------------------------------------------------CC subroutine getcm2p(iat,iatk,jat,jatk,sp2) implicit double precision(a-h,o-z) #include "divcon.dim" #include "divcon.h" dimension sp2(maxatm) i = min(iatk,jatk) j = max(iatk,jatk) c = 0.0 d = 0.0 if (i.eq.1) then C . H atom if (j.eq.6) then C . overlap with C atom c = 0.003 elseif (j.eq.7) then C . overlap with N atom c = 0.274 elseif (j.eq.8) then C . overlap with O atom c = 0.185 elseif (j.eq.14) then C . overlap with Si atom c = -0.021 elseif (j.eq.15) then C . overlap with P atom d = 0.253 elseif (j.eq.16) then C . overlap with S atom c = 0.089 endif elseif (i.eq.6) then C . C atom if (j.eq.7) then C . overlap with N atom c = -0.022 d = 0.156 elseif (j.eq.8) then C . overlap with O atom c = 0.025 d = 0.016 elseif (j.eq.9) then C . overlap with F atom d = 0.025 elseif (j.eq.14) then C . overlap with Si atom c = -0.107 elseif (j.eq.15) then C . overlap with P atom d = 0.082 elseif (j.eq.16) then C . overlap with S atom c = -0.033 d = 0.112 elseif (j.eq.17) then C . overlap with Cl atom d = 0.117 elseif (j.eq.35) then C . overlap with Br atom d = 0.040 elseif (j.eq.53) then C . overlap with I atom d = -0.032 endif elseif (i.eq.7) then C . N atom if (j.eq.8) then C . overlap with O atom c = -0.030 d = -0.043 endif elseif (i.eq.8) then C . O atom if (j.eq.15) then C . overlap with P atom d = 0.181 elseif (j.eq.16) then C . overlap with S atom d = 0.056 endif elseif (i.eq.9) then C . F atom if (j.eq.15) then C . overlap with P atom d = 0.244 endif elseif (i.eq.15) then C . P atom if (j.eq.16) then C . overlap with S atom d = -0.087 endif endif t = sp2(jat)*(d + c*sp2(jat)) if (i.ne.iatk) t = -t atchg3(iat) = atchg3(iat) + t end