c************************************************************** c************************************************************** c Segments from: c c MOLECULAR MECHANICS SUBPROGRAMS FOR MOLECULES c copyright (c) 1989 by Tamar Schlick c c************************************************************** c************************** SUBROUTINES *********************** subroutine vcross(a,b,cross) c vector version subroutine for a cross product implicit double precision(a-h,o-z) dimension a(3),b(3),cross(3) cross(1)=a(2)*b(3)-a(3)*b(2) cross(2)=a(3)*b(1)-b(3)*a(1) cross(3)=a(1)*b(2)-a(2)*b(1) return end c ******************************************************************* subroutine cosba( i, j, ii, k,l, m,nn, nout, co,dco1,dco2,ddco) c c Compute COS(THETA), where THETA is a Bond Angle, and its derivatives c c INPUT I,J,II - Triplet of atom indices making up THETA c ----- K,L,M,NN - Four coordinate indices (see below) c NOUT - Integer to specify desired ouput: c 0 - CO = COS(THETA) only c 1 - CO and DCO1 - its 1st partial der. wrt XX(K,L) c 2 - CO, its two 1st partial der. wrt XX(K,L) and c XX(M,NN), (DCO1 and DCO2 resp.), and DDCO - c the assoc. 2nd partial der. c c OUTPUT CO - COS(THETA) c ------ DCO1,DCO2 - 1st partial derivatives of CO c DDCO - 2nd partial derivative of CO c implicit double precision(a-h,o-z) parameter(maxa=200) character iatom(maxa)*4 common/coord/xx(maxa,3),iatom,id(maxa) dimension a(3), b(3), it(2), in(2) dimension dera(2),derb(2),daa(2),dbb(2),dnum(2),dden(2),dcosa(2) integer q double precision num, den double precision dot dot(z1,z2,z3,y1,y2,y3)=z1*y1 + z2*y2 +z3*y3 c c external function used: DER1(i,j,k) and DELTA(i,j) c DELTA(i,j) = 1. if i=j, 0. else c DER1(i,j,k) = 1. if k=i, -1. if k=j c do 10 ind = 1, 3 a(ind) = xx( i,ind) - xx( j,ind) b(ind) = xx(ii,ind) - xx( j,ind) 10 continue ab = dot ( a(1), a(2), a(3), b(1), b(2), b(3) ) aa = dot ( a(1), a(2), a(3), a(1), a(2), a(3) ) bb = dot ( b(1), b(2), b(3), b(1), b(2), b(3) ) num = ab den = sqrt ( aa * bb ) co = num / den if (nout.eq.0) return c COMPUTE derivatives of the NUMerator and DENominator c Derivatives of DEN are computed from derivatives of aa and bb: c d1 ( Y ** 1/2 ) = 1/2 (Y ** -1/2) d1(Y) = d1(Y) / (2. * ( Y **1/2)) c where Y = aa * bb, Y ** 1/2 = DEN. dY is calculated as c derivative of a product. Second derivatives of DEN are computed c by differentiating again this relation. This produces the formula: c d1d2 ( DEN ) = d1d2(Y) / (2*DEN) - ( d1(DEN)*d2(DEN) ) / DEN. c d1d2(Y) has 4 terms since it comes from the sum of two product terms. dco1 = 0.0 match1 = 0 match2 = 0 if ( (k.eq.i).or.(k.eq.j).or.(k.eq.ii) ) match1 = 1 if ( (m.eq.i).or.(m.eq.j).or.(m.eq.ii) ) match2 = 1 if ((nout.eq.1).and. (match1.eq.0) ) return dco2 = 0.0 ddco = 0.0 if ((nout.eq.2).and. (match1+match2.eq.0) ) return it(1) = k in(1) = l it(2) = m in(2) = nn do 20 q = 1, nout dera(q) = der1 (i , j, it(q)) derb(q) = der1 (ii, j, it(q)) dnum(q) = dera(q) * b(in(q)) + a(in(q)) * derb(q) daa(q) = 2. * dera(q) * a(in(q)) dbb(q) = 2. * derb(q) * b(in(q)) dden(q) = ( aa * dbb(q) + daa(q) * bb ) / (2. * den) dcosa(q) = (dnum(q) - (dden(q) * co) ) / den 20 continue dco1 = dcosa(1) if (nout.eq.1) return dco2 = dcosa(2) del = delta (nn,l) ddnum = del * (dera(1) * derb(2) + dera(2)* derb(1) ) ddaa = del * ( 2. * dera(1) * dera(2) ) ddbb = del * ( 2. * derb(1) * derb(2) ) ddden = ( (aa*ddbb + bb*ddaa + daa(1)*dbb(2) + daa(2)*dbb(1)) # / (2. * den) - ( dden(1) * dden(2)) / (den) ) ddco = ( ddnum - (ddden*co) - (dden(1)*dcosa(2)) # - (dden(2)*dcosa(1)) ) / den c NOUT = 2, return all 4 values return end c*********************************************************************** subroutine cosda( i,j, ii,jj, k,l,m,nn, nout, co,dco1,dco2,ddco) c c Compute COS(TAU), where TAU is a Dihedral Angle, and its derivatives c c INPUT I,J,II,JJ - Quadruplet of atom indices making up TAU. c ----- K,L,M,NN - Four coordinate indices (see below) c NOUT - Integer to specify desired ouput: c 0 - CO = COS(TAU) only c 1 - CO and DCO1 - its 1st partial der. wrt XX(K,L) c 2 - CO, its two 1st partial der. wrt XX(K,L) and c XX(M,NN), (DCO1 and DCO2 resp.), and DDCO - c the assoc. 2nd partial der. c c OUTPUT CO - COS(TAU) c ------ DCO1,DCO2 - 1st partial derivatives of CO c DDCO - 2nd partial derivative of CO c implicit double precision(a-h,o-z) parameter(maxa=200) character iatom(maxa)*4 common/coord/xx(maxa,3),iatom,id(maxa) dimension a(3), b(3), c(3), it(2), in(2) integer q dimension dera(2),derb(2),derc(2), dab(2),dbc(2), dac(2), daa(2), # dbb(2), dcc(2), dnum(2), daxb(2), dbxc(2), dden(2), dcosa(2) double precision num double precision dot dot(z1,z2,z3,y1,y2,y3)=z1*y1 + z2*y2 +z3*y3 c c external functions used: DER1(i,j,k), DELTA(i,j) c DELTA(i,j) = 1. if i=j, 0. else c DER1(i,j,k) = 1. if k=i, -1. if k=j c do 10 ind = 1, 3 a(ind) = xx( j,ind) - xx( i,ind) b(ind) = xx(ii,ind) - xx( j,ind) c(ind) = xx(jj,ind) - xx(ii,ind) 10 continue ab = dot ( a(1), a(2), a(3), b(1), b(2), b(3) ) bc = dot ( b(1), b(2), b(3), c(1), c(2), c(3) ) ac = dot ( a(1), a(2), a(3), c(1), c(2), c(3) ) aa = dot ( a(1), a(2), a(3), a(1), a(2), a(3) ) bb = dot ( b(1), b(2), b(3), b(1), b(2), b(3) ) cc = dot ( c(1), c(2), c(3), c(1), c(2), c(3) ) c Let "x" and "." denote the vector cross and dot operators, resp. Then, c NUMerator = (a x b ) . (b x c) = (a . b) (b . c) - (a . c ) (b . b) c DEN = { || a x b || * || b x c || } = { AXB * BXC } **1/2 where c AXB = (a x b) .( a x b) = (a . a) (b . b) - (a . b )**2, c BXC = (b x c) .( b x c) = (b . b) (c . c) - (b . c )**2, c axb = (aa * bb) - (ab * ab ) bxc = (bb * cc) - (bc * bc ) num = ( (ab * bc) - (ac * bb) ) den = sqrt (axb * bxc) co = num / den if (nout.eq.0) return c COMPUTE derivatives of the NUMerator and DENominator c Derivatives of DEN are computed from derivatives of axb and bxc: c d1 ( Y ** 1/2 ) = 1/2 (Y ** -1/2) d1(Y) = d1(Y) / (2. * ( Y **1/2)) c where Y = axb * bxc, Y ** 1/2 = DEN. Derivatives of axb and bxc are c calculated from the derivatives of the dot product factors, c and then dY is calculated as derivative of a product. c Second derivatives of DEN are computed by differentiating again the c relation above. This produces the formula: c d1d2 ( DEN ) = d1d2(Y) / (2*DEN) - ( d1(DEN)*d2(DEN) ) / DEN. c d1d2(Y) has 4 terms since it comes from the sum of two product terms. dco1 = 0.0 match1 = 0 match2 = 0 if ( (k.eq.i).or.(k.eq.j).or.(k.eq.ii).or.(k.eq.jj) ) match1 = 1 if ( (m.eq.i).or.(m.eq.j).or.(m.eq.ii).or.(m.eq.jj) ) match2 = 1 if ((nout.eq.1).and. (match1.eq.0) ) return dco2 = 0.0 ddco = 0.0 if ((nout.eq.2).and. ((match1+match2).eq.0) ) return it(1) = k in(1) = l it(2) = m in(2) = nn do 20 q = 1, nout dera(q) = der1(j, i, it(q) ) derb(q) = der1(ii,j, it(q) ) derc(q) = der1(jj,ii,it(q) ) dab(q) = dera(q) * b(in(q)) + a(in(q)) * derb(q) dbc(q) = derb(q) * c(in(q)) + b(in(q)) * derc(q) dac(q) = dera(q) * c(in(q)) + a(in(q)) * derc(q) daa(q) = 2. * dera(q) * a(in(q)) dbb(q) = 2. * derb(q) * b(in(q)) dcc(q) = 2. * derc(q) * c(in(q)) dnum(q) = (ab* dbc(q) + dab(q)* bc) - ( ac* dbb(q) + dac(q)* bb) daxb(q) = (aa* dbb(q) + daa(q) * bb) - ( 2. * ab * dab(q)) dbxc(q) = (bb* dcc(q) + dbb(q) * cc) - ( 2. * bc * dbc(q)) dden(q) = ( axb * dbxc(q) + daxb(q) * bxc ) / ( 2. * den) dcosa(q) = (dnum(q) - (dden(q) * co) ) / den 20 continue dco1 = dcosa(1) if (nout.eq.1) return dco2 = dcosa(2) del = delta (nn,l) ddaa = del * ( 2. * dera(1) * dera(2) ) ddbb = del * ( 2. * derb(1) * derb(2) ) ddcc = del * ( 2. * derc(1) * derc(2) ) ddab = del * ( dera(1) * derb(2) + dera(2) * derb(1) ) ddbc = del * ( derb(1) * derc(2) + derb(2) * derc(1) ) ddac = del * ( dera(1) * derc(2) + dera(2) * derc(1) ) ddnum = (ab * ddbc + dab(2)*dbc(1) + dab(1)*dbc(2) + ddab * bc ) # -(ac * ddbb + dac(2)*dbb(1) + dac(1)*dbb(2) + ddac * bb ) ddaxb = (aa * ddbb + daa(2)*dbb(1) + daa(1)*dbb(2) + ddaa * bb ) # -2. * (ab * ddab + dab(1) * dab(2) ) ddbxc = (bb * ddcc + dbb(2)*dcc(1) + dbb(1)*dcc(2) + ddbb * cc ) # -2. * (bc * ddbc + dbc(1) * dbc(2) ) ddden=((axb*ddbxc + daxb(1)*dbxc(2) + daxb(2)*dbxc(1) + ddaxb*bxc) # / (2. * den) - ( dden(1) * dden(2)) / (den) ) ddco = ( ddnum - (ddden*co) - (dden(1)*dcosa(2)) # - (dden(2)*dcosa(1)) ) / den c NOUT = 2, return all 4 values return end c ********************************************* c ***************** FUNCTIONS ***************** double precision function der1(i,j,k) c der1 is the partial derivative of xx(i,l) - xx(j,l) c wrt coordinate xx(k,l) implicit double precision(a-h,o-z) der1=0.0d0 if (k.eq.i) der1= 1.0d0 if (k.eq.j) der1= -1.0d0 return end c********************************************** double precision function delta(i,j) c the usual DELTA function delta = 0.0 if (i.eq.j) delta = 1.0 return end c********************************************** double precision function length(x,y,z) c length of a vector of components x,y,z implicit double precision(a-h,o-z) length=sqrt(x**2+y**2+z**2) return end c********************************************** double precision function vlen(a) c length of a vector a implicit double precision(a-h,o-z) dimension a(3) vlen=sqrt(a(1)**2+a(2)**2+a(3)**2) return end c********************************************** double precision function dtermt(a,b,c) c determinant of the 3*3 matrix with ROW vectors a,b,c implicit double precision(a-h,o-z) double precision a(3),b(3),c(3) dtermt=a(1)* (b(2)*c(3) - c(2)*b(3)) # -a(2)* (b(1)*c(3) - c(1)*b(3)) # +a(3)* (b(1)*c(2) - c(1)*b(2)) return end c********************************************** double precision function ssign(x) c sign of a quantity implicit double precision(a-h,o-z) if (x.gt.0.0) ssign=1.0 if (x.eq.0.0) ssign=0.0 if (x.lt.0.0) ssign=-1.0 return end c**********************************************