[CPMD-list] Torsion constraints
Ari P Seitsonen
Ari.P.Seitsonen at iki.fi
Tue Apr 27 08:08:32 CEST 2004
Dear Erik,
Here's a version of the routines (in 'constr.F') which _should_ now be
correct, I'm not 100 % sure that this file alone is enough to fix the bug
but I think so.
Greetings,
apsi
-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-
Ari Paavo Seitsonen / Ari.P.Seitsonen at iki.fi / http://www.iki.fi/~apsi/
Tel +41 1 635 44 97 / Fax +41 1 635 68 38 / GSM +41 79 719 09 35
Anschrift: Physikalisch Chemisches Institut (PCI), Universität Zürich (UniZh)
Indirizzo: Winterthurerstraße 190, CH-8057 Zürich
Address: Schweiz / Svizzera / Suisse / Svizra / Switzerland
On Mon, 26 Apr 2004, Erik Santiso wrote:
>
> Good! At least I know I wasn't defining the constraints wrong :)
>
> I've been able to get it to work by setting the torsion angle to some
> value close to zero but not zero (say 0.01). Sometimes it still jumps and
> shows a huge value of GNMAX, but at the end it converges. Guess I'll have
> to check the numbers again when the new version comes out.
>
> Thanks for your help!
>
> Erik.
>
> -------------------------------------
> Vir prudens non contra ventum mingit.
>
> >
> > Dear Erik,
> >
> > There's indeed a bug in the CPMD version 3.7.*, and it's exactly the
> > phase of the angle (-1 degree becomes 359 or vice versa). It should be
> > fixed in the new version (out any day now...!! I hope). Sorry for the
> > inconvenience.
> >
> > GNMAX is the largest simple component of ionic forces, however
> > what's printed with the coordinates are the forces without/before the
> > applying the constraints, that's why GNMAX is huge even though the forces
> > in the output seem to be small.
> >
> > Greetings,
> >
> > apsi
> >
> > -=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-
> > Ari Paavo Seitsonen / Ari.P.Seitsonen at iki.fi / http://www.iki.fi/~apsi/
> > Tel +41 1 635 44 97 / Fax +41 1 635 68 38 / GSM +41 79 719 09 35
> > Anschrift: Physikalisch Chemisches Institut (PCI), Universität Zürich
> > (UniZh)
> > Indirizzo: Winterthurerstraße 190, CH-8057 Zürich
> > Address: Schweiz / Svizzera / Suisse / Svizra / Switzerland
> >
> > On Sun, 18 Apr 2004, Erik Santiso wrote:
> >
> >> Hi,
> >>
> >> I've been doing some geometry optimizations with torsion constraints,
> >> and
> >> for most cases thinks seem to work well. However, if I set a torsion
> >> angle
> >> to be 0 degrees (or 360), the geometry optimization starts doing funny
> >> things. Here's an example. The constraint in this case is that the
> >> torsion
> >> between atoms 7, 8, 9, 10 should be 0.0:
> >>
> >> ATOM COORDINATES GRADIENTS (-FORCES)
> >> 1 H 0.6369 11.2421 7.3588 0.132E-03 -0.197E-04 0.250E-03
> >> 2 H 2.7954 13.9975 7.3282 -0.144E-03 0.423E-04 -0.288E-03
> >> 3 H -3.0319 13.9637 7.3442 0.110E-03 0.155E-04 -0.891E-04
> >> 4 H -3.0301 18.4817 7.3516 -0.723E-04 0.549E-04 0.348E-03
> >> 5 H 0.6320 21.2032 7.3462 0.238E-04 -0.108E-03 -0.193E-03
> >> 6 H 2.7974 18.4518 7.3616 -0.199E-04 -0.110E-03 0.265E-03
> >> 7 C 0.8615 13.2867 7.3435 0.171E-03 -0.140E-03 0.750E-03
> >> 8 C -1.1573 14.8298 7.3389 -0.341E-03 0.187E-03 -0.477E-03
> >> 9 C -1.1552 17.6146 7.3452 0.332E-03 -0.267E-03 0.557E-03
> >> 10 C 0.8619 19.1595 7.3501 -0.956E-05 0.325E-03 0.346E-03
> >> ****************************************************************
> >> *** TOTAL STEP NR. 166 GEOMETRY STEP NR. 24 ***
> >> *** GNMAX= 1.696295E-02 [8.02E-03] ETOT= -26.304652 ***
> >> *** GNORM= 4.848471E-03 DETOT= -1.048E-05 ***
> >> *** CNSTR= 1.745681E-02 TCPU= 31.09 ***
> >> ****************************************************************
> >> 1 2.674E-04 4.708E-05 -26.304543 1.082E-04 4.04
> >> 2 1.048E-04 1.539E-05 -26.304644 -1.002E-04 4.07
> >> 3 5.823E-05 8.073E-06 -26.304651 -7.269E-06 4.05
> >> 4 5.257E-05 4.483E-06 -26.304654 -3.038E-06 4.10
> >> 5 2.776E-05 2.903E-06 -26.304654 -6.502E-07 4.10
> >> 6 1.711E-05 1.621E-06 -26.304655 -4.005E-07 4.90
> >> 7 1.186E-05 6.911E-07 -26.304655 -2.197E-07 4.78
> >> 8 6.332E-06 3.436E-07 -26.304655 -6.332E-09 4.85
> >>
> >> RESTART INFORMATION WRITTEN ON FILE ./RESTART.1
> >>
> >> ATOM COORDINATES GRADIENTS (-FORCES)
> >> 1 H 0.6355 11.2413 7.3588 0.526E-04 -0.877E-04 0.409E-03
> >> 2 H 2.7954 13.9970 7.3292 -0.168E-03 0.962E-04 -0.229E-04
> >> 3 H -3.0314 13.9634 7.3479 0.132E-03 -0.633E-05 -0.574E-04
> >> 4 H -3.0297 18.4816 7.3509 -0.160E-03 0.100E-03 0.300E-03
> >> 5 H 0.6309 21.2036 7.3458 0.320E-04 -0.158E-03 -0.187E-03
> >> 6 H 2.7984 18.4542 7.3634 0.473E-04 -0.805E-04 0.287E-03
> >> 7 C 0.8617 13.2859 7.3373 0.252E-03 -0.145E-03 0.562E-04
> >> 8 C -1.1569 14.8295 7.3402 -0.414E-03 0.118E-03 -0.197E-03
> >> 9 C -1.1545 17.6146 7.3461 0.506E-03 -0.255E-03 0.583E-03
> >> 10 C 0.8621 19.1602 7.3504 -0.729E-04 0.368E-03 0.278E-03
> >> ****************************************************************
> >> *** TOTAL STEP NR. 174 GEOMETRY STEP NR. 25 ***
> >> *** GNMAX= 3.285935E+01 [6.14E-03] ETOT= -26.304655 ***
> >> *** GNORM= 9.384478E+00 DETOT= -3.556E-06 ***
> >> *** CNSTR= 3.285993E+01 TCPU= 34.97 ***
> >> ****************************************************************
> >> 1 2.379E-04 4.148E-05 -26.304568 8.690E-05 4.05
> >> 2 9.553E-05 1.362E-05 -26.304642 -7.339E-05 4.02
> >> 3 5.412E-05 7.187E-06 -26.304650 -8.253E-06 4.04
> >> 4 4.862E-05 3.975E-06 -26.304651 -1.259E-06 4.21
> >> 5 2.584E-05 2.507E-06 -26.304652 -7.716E-07 4.10
> >> 6 1.562E-05 1.417E-06 -26.304652 -3.099E-07 4.83
> >> 7 1.029E-05 5.705E-07 -26.304652 -7.033E-08 4.83
> >> 8 5.174E-06 2.462E-07 -26.304652 -3.037E-08 4.79
> >>
> >> RESTART INFORMATION WRITTEN ON FILE ./RESTART.1
> >>
> >> ATOM COORDINATES GRADIENTS (-FORCES)
> >> 1 H 0.6367 11.2421 7.3587 0.158E-03 0.159E-04 0.284E-03
> >> 2 H 2.7959 13.9965 7.3279 -0.479E-04 0.339E-04 -0.247E-03
> >> 3 H -3.0317 13.9639 7.3445 0.308E-04 0.242E-04 -0.100E-03
> >> 4 H -3.0299 18.4821 7.3517 0.805E-04 0.338E-04 0.350E-03
> >> 5 H 0.6315 21.2030 7.3463 0.164E-04 -0.373E-04 -0.187E-03
> >> 6 H 2.7977 18.4529 7.3624 -0.468E-05 -0.183E-04 0.273E-03
> >> 7 C 0.8614 13.2866 7.3426 -0.168E-03 -0.316E-04 0.715E-03
> >> 8 C -1.1568 14.8297 7.3389 0.467E-04 0.347E-04 -0.483E-03
> >> 9 C -1.1557 17.6147 7.3453 0.256E-04 -0.184E-03 0.562E-03
> >> 10 C 0.8618 19.1591 7.3502 0.529E-04 0.812E-04 0.308E-03
> >> ****************************************************************
> >> *** TOTAL STEP NR. 182 GEOMETRY STEP NR. 26 ***
> >> *** GNMAX= 1.532222E-02 [5.30E-03] ETOT= -26.304652 ***
> >> *** GNORM= 4.377040E-03 DETOT= 2.814E-06 ***
> >> *** CNSTR= 1.581241E-02 TCPU= 34.95 ***
> >> ****************************************************************
> >> 1 8.455E-06 2.055E-06 -26.304652 2.594E-07 4.12
> >> 2 6.064E-06 6.952E-07 -26.304652 -2.099E-07 4.85
> >>
> >> RESTART INFORMATION WRITTEN ON FILE ./RESTART.1
> >>
> >> In the middle step (25), it appears that the constraint has a difference
> >> of 32.8, but the torsion angle is actually much smaller (the coordinates
> >> of the 4 atoms are pretty much equal to those in the other two steps).
> >> Any
> >> idea what is causing this? One thing that comes to mind is the
> >> discontinuity of the torsion angle at 0 (360), could this be the cause
> >> (and if it is, is there any way around it?).
> >>
> >> I also have another question. What are exactly the quantities in the
> >> GNMAX
> >> line of the output at each geometry step? Before I thought it was the
> >> maximum nuclear gradient, but in this case it is clearly not.
> >>
> >> Thanks a lot!
> >>
> >> Erik.
> >>
> >> -------------------------------------
> >> Vir prudens non contra ventum mingit.
> >> _______________________________________________
> >> CPMD-list mailing list
> >> CPMD-list at cpmd.org
> >> http://www.cpmd.org/mailman/listinfo/cpmd-list
> >>
> >
>
>
>
> _______________________________________________
> CPMD-list mailing list
> CPMD-list at cpmd.org
> http://www.cpmd.org/mailman/listinfo/cpmd-list
>
-------------- next part --------------
C----------------STRETCH-----------------------------------------------C
subroutine funcr(fr,r,R0,x,y)
c FUNCTION: R*R - R0*R0
implicit none
c Argument variables
real*8 fr, r, r0, x(3), y(3)
c Local variables
real*8 t11, t3, t6, t9
t3 = (x(1)-y(1))**2
t6 = (x(2)-y(2))**2
t9 = (x(3)-y(3))**2
t11 = t3+t6+t9
fr = t11-R0*R0
r = sqrt(t11)
return
end
C----------------DISTANCE----------------------------------------------C
subroutine funcd(fd,d,R0,x,y)
c FUNCTION: R - R0
implicit none
c Argument variables
real*8 d, fd, r0, x(3), y(3)
c Local variables
real*8 t3, t6, t9
t3 = (x(1)-y(1))**2
t6 = (x(2)-y(2))**2
t9 = (x(3)-y(3))**2
d = sqrt(t3+t6+t9)
fd = d-R0
return
end
C----------------DIFFERENCE OF DISTANCES-------------------------------C
subroutine funcdd(fd,d,R0,x,y,z)
c FUNCTION: dR - R0
implicit none
c Argument variables
real*8 d, fd, r0, x(3), y(3), z(3)
c Local variables
real*8 d1, d2, t1, t2, t3, t4, t5, t6
t1 = (x(1)-y(1))**2
t2 = (x(2)-y(2))**2
t3 = (x(3)-y(3))**2
t4 = (y(1)-z(1))**2
t5 = (y(2)-z(2))**2
t6 = (y(3)-z(3))**2
d1 = sqrt(t1+t2+t3)
d2 = sqrt(t4+t5+t6)
d = d1-d2
fd = d-R0
return
end
C----------------ANGLE-------------------------------------------------C
subroutine funct(ft,t,T0,x,y,z)
c FUNCTION: COS(theta) - COS(theta0)
implicit none
c Argument variables
real*8 ft, t, t0, x(3), y(3), z(3)
c Local variables
real*8 epsilon
parameter (epsilon=1.D-12)
real*8 ac, caa, cab, cbb, t1, t2, t3, t4, t5, t6
t1 = (x(1)-y(1))
t2 = (x(2)-y(2))
t3 = (x(3)-y(3))
t4 = (y(1)-z(1))
t5 = (y(2)-z(2))
t6 = (y(3)-z(3))
caa = t1*t1+t2*t2+t3*t3
cab = t1*t4+t2*t5+t3*t6
cbb = t4*t4+t5*t5+t6*t6
if(caa.lt.epsilon.or.cbb.lt.epsilon) then
C 2 points at the same place (T.D. 16/12/1999).
C We give the desired value (T0) for satisfied constraint.
C ac=-cos(T0)
ft=0.D0
t = T0
else
ac = cab/sqrt(caa*cbb)
if(ac.lt.-1.D0) ac=-1.D0
if(ac.gt.1.D0) ac=1.D0
ft = ac+cos(T0)
t = acos(-ac)
endif
return
end
C----------------DIHEDRAL----------------------------------------------C
subroutine funco(fo,o,O0,x,y,z,w,sign0)
c FUNCTION: phi - phi0
implicit none
c Argument variables
real*8 fo, o, o0, w(3), x(3), y(3), z(3),sign0
c Local variables
real*8 epsilon
parameter (epsilon=1.D-12)
real*8 ac, caa, cab, cac, cbb, cbc, ccc, dab, dbc, t1, t2, t3,
& t4, t5, t6, t7, t8, t9
real*8 a(3),b(3),c(3),avb(3),bvc(3),vv(3),s1,pi
integer i
pi=acos(-1.)
do i=1,3
a(i) = (x(i)-y(i))
b(i) = (z(i)-y(i))
enddo
call normalize(a)
call normalize(b)
call vecprod(a,b,avb)
call normalize(avb)
do i=1,3
b(i) = (y(i)-z(i))
c(i) = (w(i)-z(i))
enddo
call normalize(b)
call normalize(c)
call vecprod(b,c,bvc)
call normalize(bvc)
call getscal(avb,bvc,ac)
call vecprod(avb,bvc,vv)
do i=1,3
b(i) = (z(i)-y(i))
enddo
call getscal(b,vv,s1)
o=acos(ac)
! a = x-y
! b = z-y (or y-z)
! c = w-z
! avb = a x -b / || a x -b ||
! bvc = b x c / || b x c ||
! ac = avb . bvc = ( a x -b ) . ( b x c ) / ( || a x -b || * || b x c || )
! = - ( ( a x b ) x b ) . c / ( || a x -b || * || b x c || )
! = [ ( b . b ) ( a . c ) - ( a . b ) ( b . c ) ] / ( || || * || || )
sign0=+1.
if(s1.lt.0.d0)then
sign0=-1.
o=-o
endif
fo=mod(20000.0*pi+o-o0,2.0*pi)
if ( fo .gt. +pi ) fo = fo - 2.*pi
o=o0+fo
return
end
!
subroutine normalize(b)
real*8 b(3),norm
norm=sqrt(b(1)**2+b(2)**2+b(3)**2)
b(1)=b(1)/norm
b(2)=b(2)/norm
b(3)=b(3)/norm
return
end
C----------------OUT OF PLANE------------------------------------------C
SUBROUTINE funcp(diff,teta,teta0,x,y,z,w,d)
c FUNCTION: phi - phi0
implicit none
c Argument variables
real*8 diff,teta,teta0,w(3),x(3),y(3),z(3),d(12)
c Local variables
real*8 epsilon
parameter (epsilon=1.D-12)
real*8 sinus,cosine,t(3),q(3),nn(3),o(3),pi,sign0,t_norm
real*8 ovn(3),num,den,const,o_norm,nn_norm,nn2,o2,scal
integer i
sign0=+1.
pi=acos(-1.D0)
do i=1,3
t(i) = (x(i)-y(i))
q(i) = (z(i)-y(i))
enddo
t_norm = t(1)*t(1)+t(2)*t(2)+t(3)*t(3)
t_norm = dsqrt(t_norm)
call vecprod(t,q,nn)
nn2 = nn(1)*nn(1)+nn(2)*nn(2)+nn(3)*nn(3)
nn_norm = dsqrt(nn2)
do i=1,3
o(i) = (w(i)-y(i))
enddo
o2 = o(1)*o(1)+o(2)*o(2)+o(3)*o(3)
o_norm = dsqrt(o2)
call getscal(o,t,scal)
call getscal(o,nn,cosine)
num = cosine
den = 1.0D0/(o_norm*nn_norm)
cosine = cosine*den
sinus = dsqrt(1-cosine*cosine)
teta = acos(cosine)
if(scal .lt. 0.d0) then
sign0=-1.d0
teta=-teta+2.d0*pi
endif
diff=mod(20000*pi+teta0,2.*pi)-teta
const = -1.0D0/sinus*sign0
d(1) = const*den*(-q(3)*o(2)+q(2)*o(3)
& -num*(-q(3)*nn(2)+q(2)*nn(3))/nn2)
d(2) = const*den*(-q(1)*o(3)+q(3)*o(1)
& -num*(-q(1)*nn(3)+q(3)*nn(1))/nn2)
d(3) = const*den*(-q(2)*o(1)+q(1)*o(2)
& -num*(-q(2)*nn(1)+q(1)*nn(2))/nn2)
d(4) = const*den*(-nn(1)-t(3)*o(2)+q(3)*o(2)-q(2)*o(3)+t(2)*o(3)
& -num*(-o(1)/o2+
& ((-t(3)+q(3))*nn(2)+(-q(2)+t(2))*nn(3))/nn2))
d(5) = const*den*(-nn(2)-q(3)*o(1)+t(3)*o(1)-t(1)*o(3)+q(1)*o(3)
& -num*(-o(2)/o2+
& ((-q(3)+t(3))*nn(1)+(-t(1)+q(1))*nn(3))/nn2))
d(6) = const*den*(-nn(3)-t(2)*o(1)+q(2)*o(1)-q(1)*o(2)+t(1)*o(2)
& -num*(-o(3)/o2+
& ((-t(2)+q(2))*nn(1)+(-q(1)+t(1))*nn(2))/nn2))
d(7) = const*den*(t(3)*o(2)-t(2)*o(3)
& -num*(t(3)*nn(2)-t(2)*nn(3))/nn2)
d(8) = const*den*(t(1)*o(3)-t(3)*o(1)
& -num*(t(1)*nn(3)-t(3)*nn(1))/nn2)
d(9) = const*den*(t(2)*o(1)-t(1)*o(2)
& -num*(t(2)*nn(1)-t(1)*nn(2))/nn2)
d(10) = const*den*(nn(1)-num*o(1)/o2)
d(11) = const*den*(nn(2)-num*o(2)/o2)
d(12) = const*den*(nn(3)-num*o(3)/o2)
return
end
C----------------------------------------------------------------------C
subroutine diffd(dR,x,y)
implicit none
c Argument variables
real*8 dR(6), x(3),y(3)
c Local variables
real*8 epsilon
parameter (epsilon=1.D-12)
real*8 r, t3, t6, t9
t3 = (x(1)-y(1))
t6 = (x(2)-y(2))
t9 = (x(3)-y(3))
r = sqrt(t3*t3+t6*t6+t9*t9)
IF(r.lt.epsilon) then
dR(1) = 0.D0
dR(2) = 0.D0
dR(3) = 0.D0
dR(4) = 0.D0
dR(5) = 0.D0
dR(6) = 0.D0
else
dR(1) = t3/r
dR(2) = t6/r
dR(3) = t9/r
dR(4) = -t3/r
dR(5) = -t6/r
dR(6) = -t9/r
endif
return
end
C----------------------------------------------------------------------C
subroutine diffdd(dR,x,y,z)
implicit none
c Argument variables
real*8 dR(9), x(3), y(3), z(3)
c Local variables
real*8 epsilon
parameter (epsilon=1.D-12)
real*8 r1, r2, t1, t2, t3, t4, t5, t6
t1 = (x(1)-y(1))
t2 = (x(2)-y(2))
t3 = (x(3)-y(3))
t4 = (y(1)-z(1))
t5 = (y(2)-z(2))
t6 = (y(3)-z(3))
r1 = sqrt(t1*t1+t2*t2+t3*t3)
r2 = sqrt(t4*t4+t5*t5+t6*t6)
if(r1.lt.epsilon) then
dR(1) = 0.D0
dR(2) = 0.D0
dR(3) = 0.D0
else
dR(1) = t1/r1
dR(2) = t2/r1
dR(3) = t3/r1
endif
if(r2.lt.epsilon) then
dR(7) = 0.D0
dR(8) = 0.D0
dR(9) = 0.D0
else
dR(7) = t4/r2
dR(8) = t5/r2
dR(9) = t6/r2
endif
dR(4) = -dR(1)-dR(7)
dR(5) = -dR(2)-dR(8)
dR(6) = -dR(3)-dR(9)
return
end
C----------------------------------------------------------------------C
subroutine diffr(dR,x,y)
implicit none
c Argument variables
real*8 dR(6), x(3),y(3)
c Local variables
real*8 t3, t6, t9
t3 = (x(1)-y(1))
t6 = (x(2)-y(2))
t9 = (x(3)-y(3))
dR(1) = 2.D0*t3
dR(2) = 2.D0*t6
dR(3) = 2.D0*t9
dR(4) = -2.D0*t3
dR(5) = -2.D0*t6
dR(6) = -2.D0*t9
return
end
C----------------------------------------------------------------------C
subroutine difft(dT,x,y,z)
implicit none
c Argument variables
real*8 dT(9), x(3), y(3), z(3)
c Local variables
real*8 epsilon
parameter (epsilon=1.D-12)
real*8 caa, cab, cbb, ccc, t1, t2, t3, t4, t5, t6
integer i
t1 = (x(1)-y(1))
t2 = (x(2)-y(2))
t3 = (x(3)-y(3))
t4 = (y(1)-z(1))
t5 = (y(2)-z(2))
t6 = (y(3)-z(3))
caa = t1*t1+t2*t2+t3*t3
cab = t1*t4+t2*t5+t3*t6
cbb = t4*t4+t5*t5+t6*t6
if(caa.lt.epsilon.or.cbb.lt.epsilon) then
do i=1,9
dT(i)=0.D0
enddo
else
ccc = -1.D0/sqrt(caa*cbb)
dT(1) = ccc*(cab/caa*t1-t4)
dT(2) = ccc*(cab/caa*t2-t5)
dT(3) = ccc*(cab/caa*t3-t6)
dT(4) = -ccc*(cab/caa*t1-cab/cbb*t4+t1-t4)
dT(5) = -ccc*(cab/caa*t2-cab/cbb*t5+t2-t5)
dT(6) = -ccc*(cab/caa*t3-cab/cbb*t6+t3-t6)
dT(7) = -ccc*(cab/cbb*t4-t1)
dT(8) = -ccc*(cab/cbb*t5-t2)
dT(9) = -ccc*(cab/cbb*t6-t3)
endif
return
end
C----------------------------------------------------------------------C
subroutine diffo(dO,x,y,z,w,sign0)
implicit none
c Argument variables
real*8 dO(12), w(3), x(3), y(3), z(3),sign0
c Local variables
real*8 epsilon
parameter (epsilon=1.D-12)
real*8 caa, cab, cac, cbb, cbc, ccc, dab, dbc, ddd, t1, t2, t3,
& t4, t5, t6, t7, t8, t9
integer i
real*8 cost,ac
t1 = (x(1)-y(1))
t2 = (x(2)-y(2))
t3 = (x(3)-y(3))
t4 = (y(1)-z(1))
t5 = (y(2)-z(2))
t6 = (y(3)-z(3))
t7 = (z(1)-w(1))
t8 = (z(2)-w(2))
t9 = (z(3)-w(3))
caa = t1*t1+t2*t2+t3*t3
cbb = t4*t4+t5*t5+t6*t6
ccc = t7*t7+t8*t8+t9*t9
cab = t1*t4+t2*t5+t3*t6
cac = t1*t7+t2*t8+t3*t9
cbc = t4*t7+t5*t8+t6*t9
dab = caa*cbb-cab*cab
dbc = cbb*ccc-cbc*cbc
! || a x b ||^2 = ( a x b ) . ( a x b )
! = a . ( b x ( a x b ) )
! = a . ( ( b . b ) a - ( b . a ) b
! = ( a . a ) ( b . b ) - ( a . b ) ( a . b )
! dab = || a x b ||^2
! ac = - ( ( a . b ) ( b . c ) - ( a . c ) ( b . b ) )
! / ( || a x -b || * || b x c || )
!
! D[arccos(ac)] = -1 / sqrt ( 1 - ac^2 )
! d[arccos(ac)]/dx1
! = -1 / sqrt ( 1 - ac^2 ) * d(ac)/dx1
! = 1 / sqrt ( 1 - ac^2 ) * [ b1 ( b . c ) - ( b . b ) c1 ] / ( || || || || )
! - 1 / sqrt ( 1 - ac^2 ) * [ ( a . b ) ( b . c ) - ( a . c ) ( b . b ) ]
! * d(|| a x b ||) / dx1 * || -b x c || / ( || || || || )^2
! = 1 / sqrt ( 1 - ac^2 ) * [ b1 ( b . c ) - ( b . b ) c1 ] / ( || || || || )
! - 1 / sqrt ( 1 - ac^2 ) * [ ( a . b ) ( b . c ) - ( a . c ) ( b . b ) ]
! * d(|| a x b ||^2) / dx1 / ( 2 * || a x b ||^3 * || -b x c || )
! = 1 / sqrt ( 1 - ac^2 ) * [ b1 ( b . c ) - ( b . b ) c1 ] / ( || || || || )
! - 1 / sqrt ( 1 - ac^2 ) * [ ( a . b ) ( b . c ) - ( a . c ) ( b . b ) ]
! * ( 2 * a1 * ( b . b ) - 2 * b1 * ( a . b ) )
! / ( 2 * || a x b ||^3 * || -b x c || )
! = 1 / sqrt ( 1 - ac^2 ) * [ b1 ( b . c ) - ( b . b ) c1 ] / ( || || || || )
! - 1 / sqrt ( 1 - ac^2 ) * [ ( a . b ) ( b . c ) - ( a . c ) ( b . b ) ]
! * ( a1 * ( b . b ) - b1 * ( a . b ) ) / ( || a x b ||^3 * || -b x c || )
! = ddd * [ b1 ( b . c ) - ( b . b ) c1 ]
! - ddd * [ ( a . b ) ( b . c ) - ( a . c ) ( b . b ) ]
! * ( a1 * ( b . b ) - b1 * ( a . b ) ) / || a x b ||^2
! = ddd * [ t4 * cbc - cbb * t7 ]
! - ddd * [ cab * cbc - cac * cbb ] * ( t1 * cbb - t4 * cab ) / dab
if(dab.lt.epsilon.or.dbc.lt.epsilon) then
do i=1,12
dO(i)=0.D0
enddo
else
! LAIO
ac = -(cab*cbc-cac*cbb)/sqrt(dab*dbc)
! if(abs(ac).gt.0.99999999D0) ac=0.99999999D0
if(abs(ac).gt.0.999999999999D0) ac=0.999999999999D0
cost= 1.D0/sqrt(1.D0-ac*ac)*sign0
ddd = 1.D0/sqrt(dab*dbc)*cost
dO(1) = -ddd*
& (cbc*t4-cbb*t7-(cab*cbc-cac*cbb)/dab*(cbb*t1-cab*t4))
dO(2) = -ddd*
& (cbc*t5-cbb*t8-(cab*cbc-cac*cbb)/dab*(cbb*t2-cab*t5))
dO(3) = -ddd*
& (cbc*t6-cbb*t9-(cab*cbc-cac*cbb)/dab*(cbb*t3-cab*t6))
dO(4) = -ddd*(cbc*t1-cbc*t4+cab*t7+cbb*t7-2*cac*t4-
* (cab*cbc-cac*cbb)/dbc*(ccc*t4-cbc*t7)-
* (cab*cbc-cac*cbb)/dab*(caa*t4-cbb*t1-cab*t1+cab*t4))
dO(5) = -ddd*(cbc*t2-cbc*t5+cab*t8+cbb*t8-2*cac*t5-
* (cab*cbc-cac*cbb)/dbc*(ccc*t5-cbc*t8)-
* (cab*cbc-cac*cbb)/dab*(caa*t5-cbb*t2-cab*t2+cab*t5))
dO(6) = -ddd*(cbc*t3-cbc*t6+cab*t9+cbb*t9-2*cac*t6-
* (cab*cbc-cac*cbb)/dbc*(ccc*t6-cbc*t9)-
* (cab*cbc-cac*cbb)/dab*(caa*t6-cbb*t3-cab*t3+cab*t6))
dO(7) = -ddd*(-cbc*t1+cab*t4-cab*t7-cbb*t1+2.D0*cac*t4-
* (cab*cbc-cac*cbb)/dbc*(cbb*t7-ccc*t4-cbc*t4+cbc*t7)-
* (cab*cbc-cac*cbb)/dab*(-caa*t4+cab*t1))
dO(8) = -ddd*(-cbc*t2+cab*t5-cab*t8-cbb*t2+2.D0*cac*t5-
* (cab*cbc-cac*cbb)/dbc*(cbb*t8-ccc*t5-cbc*t5+cbc*t8)-
* (cab*cbc-cac*cbb)/dab*(-caa*t5+cab*t2))
dO(9) = -ddd*(-cbc*t3+cab*t6-cab*t9-cbb*t3+2.D0*cac*t6-
* (cab*cbc-cac*cbb)/dbc*(cbb*t9-ccc*t6-cbc*t6+cbc*t9)-
* (cab*cbc-cac*cbb)/dab*(-caa*t6+cab*t3))
dO(10)= -ddd*(-cab*t4+cbb*t1-(cab*cbc-cac*cbb)/dbc*
* (-cbb*t7+cbc*t4))
dO(11)= -ddd*(-cab*t5+cbb*t2-(cab*cbc-cac*cbb)/dbc*
* (-cbb*t8+cbc*t5))
dO(12)= -ddd*(-cab*t6+cbb*t3-(cab*cbc-cac*cbb)/dbc*
* (-cbb*t9+cbc*t6))
endif
return
end
C----------------------------------------------------------------------C
subroutine diffp(dP,x,y,z,w)
implicit none
c Argument variables
real*8 dP(12), w(3), x(3), y(3), z(3)
c Local variables
real*8 epsilon
parameter (epsilon=1.D-12)
real*8 cab, ccc, dab, ddd, t1, t2, t3, t4, t5, t6, t7, t8, t9
integer i
t1 = (x(1)-y(1))
t2 = (x(2)-y(2))
t3 = (x(3)-y(3))
t4 = (y(1)-z(1))
t5 = (y(2)-z(2))
t6 = (y(3)-z(3))
t7 = (w(1)-y(1))
t8 = (w(2)-y(2))
t9 = (w(3)-y(3))
ccc = t7*t7+t8*t8+t9*t9
cab = t7*(t2*t6-t3*t5)+t8*(t3*t4-t1*t6)+t9*(t1*t5-t2*t4)
dab = (t2*t6-t3*t5)*(t2*t6-t3*t5)+(t3*t4-t1*t6)*(t3*t4-t1*t6)+
* (t1*t5-t2*t4)*(t1*t5-t2*t4)
if(ccc.lt.epsilon.or.dab.lt.epsilon) then
do i=1,12
dP(i)=0.D0
enddo
else
ddd = 1.D0/sqrt(ccc*dab)
dP(1) = ddd*(-t8*t6+t9*t5-
* cab*(t1*(t5*t5+t6*t6)-t4*(t2*t5+t3*t6))/dab)
dP(2) = ddd*( t7*t6-t9*t4-
* cab*(t2*(t6*t6+t4*t4)-t5*(t3*t6+t1*t4))/dab)
dP(3) = ddd*(-t7*t5+t8*t4-
* cab*(t3*(t4*t4+t5*t5)-t6*(t1*t4+t2*t5))/dab)
dP(7) = -ddd*( t8*t3-t9*t2-
* cab*(t4*(t2*t2+t3*t3)-t1*(t2*t5+t3*t6))/dab)
dP(8) = -ddd*(-t7*t3+t9*t1-
* cab*(t5*(t3*t3+t1*t1)-t2*(t3*t6+t1*t4))/dab)
dP(9) = -ddd*( t7*t2-t8*t1-
* cab*(t6*(t1*t1+t2*t2)-t3*(t1*t4+t2*t5))/dab)
dP(10)= ddd*(t2*t6-t3*t5-cab*t7/ccc)
dP(11)= ddd*(t3*t4-t1*t6-cab*t8/ccc)
dP(12)= ddd*(t1*t5-t2*t4-cab*t9/ccc)
dP(4) = -dP(1)-dP(7)-dP(10)
dP(5) = -dP(2)-dP(8)-dP(11)
dP(6) = -dP(3)-dP(9)-dP(12)
endif
return
end
C----------------------------------------------------------------------C
C ==================================================================
SUBROUTINE getscal(vec1,vec2,scal)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION vec1(3), vec2(3)
C ==--------------------------------------------------------------==
scal = vec1(1)*vec2(1) + vec1(2)*vec2(2) + vec1(3)*vec2(3)
C ==--------------------------------------------------------------==
RETURN
END
C ==================================================================
C ==================================================================
SUBROUTINE getnorm(vec,dnorm)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION vec(3)
C ==--------------------------------------------------------------==
dnorm = sqrt ( vec(1)*vec(1) + vec(2)*vec(2) + vec(3)*vec(3) )
C ==--------------------------------------------------------------==
RETURN
END
C ==================================================================
C ==================================================================
SUBROUTINE normvec(vec_old,vec_norm)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION vec_old(3), vec_norm(3)
C ==--------------------------------------------------------------==
CALL getnorm(vec_old, dnorm)
vec_norm(1) = vec_old(1) / dnorm
vec_norm(2) = vec_old(2) / dnorm
vec_norm(3) = vec_old(3) / dnorm
C ==--------------------------------------------------------------==
RETURN
END
C ==================================================================
C ==================================================================
SUBROUTINE vecmul(a,vec1,vec2)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION vec1(3), vec2(3)
C ==--------------------------------------------------------------==
vec2(1) = a * vec1(1)
vec2(2) = a * vec1(2)
vec2(3) = a * vec1(3)
C ==--------------------------------------------------------------==
RETURN
END
C ==================================================================
C ==================================================================
SUBROUTINE addvec(a,vec1,b,vec2,vec3)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION vec1(3), vec2(3), vec3(3)
C ==--------------------------------------------------------------==
vec3(1)=a*vec1(1)+b*vec2(1)
vec3(2)=a*vec1(2)+b*vec2(2)
vec3(3)=a*vec1(3)+b*vec2(3)
C ==--------------------------------------------------------------==
RETURN
END
C ==================================================================
C
C ==================================================================
SUBROUTINE copyvec(vec1,vec2)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION vec1(3), vec2(3)
C ==--------------------------------------------------------------==
vec2(1) = vec1(1)
vec2(2) = vec1(2)
vec2(3) = vec1(3)
C ==--------------------------------------------------------------==
RETURN
END
C ==================================================================
C ==================================================================
SUBROUTINE zerovec(vec1)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION vec1(3)
C ==--------------------------------------------------------------==
vec1(1) = 0.0D0
vec1(2) = 0.0D0
vec1(3) = 0.0D0
C ==--------------------------------------------------------------==
C ==--------------------------------------------------------------==
RETURN
END
C ==================================================================
C ==================================================================
SUBROUTINE vecprod(vec1,vec2,vec3)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION vec1(3), vec2(3), vec3(3)
C ==--------------------------------------------------------------==
vec3(1) = vec1(2) * vec2(3) - vec1(3) * vec2(2)
vec3(2) = vec1(3) * vec2(1) - vec1(1) * vec2(3)
vec3(3) = vec1(1) * vec2(2) - vec1(2) * vec2(1)
C ==--------------------------------------------------------------==
RETURN
END
C ==================================================================
C ==================================================================
C ==================================================================
SUBROUTINE vecrotz(angle,vec1,vec2)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION vec1(3), vec2(3)
REAL*8 R(3,3)
C ==--------------------------------------------------------------==
DO i=1,3
DO j=1,3
R(i,j)=0.0D0
END DO
vec2(i)=0.0D0
R(i,i)=1.0D0
END DO
R(1,1)=+cos(angle)
R(2,2)=+cos(angle)
R(1,2)=-sin(angle)
R(2,1)=+sin(angle)
DO i=1,3
DO j=1,3
vec2(I)=vec2(I)+R(I,J)*vec1(J)
END DO
END DO
C ==--------------------------------------------------------------==
RETURN
END
C ==================================================================
C
C ==================================================================
SUBROUTINE vecroty(angle,vec1,vec2)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION vec1(3), vec2(3)
REAL*8 R(3,3)
C ==--------------------------------------------------------------==
DO i=1,3
DO j=1,3
R(i,j)=0.0D0
END DO
vec2(i)=0.0D0
R(i,i)=1.0D0
END DO
R(1,1)=+cos(angle)
R(3,3)=+cos(angle)
R(1,3)=+sin(angle)
R(3,1)=-sin(angle)
DO i=1,3
DO j=1,3
vec2(I)=vec2(I)+R(I,J)*vec1(J)
END DO
END DO
C ==--------------------------------------------------------------==
RETURN
END
C ==================================================================
C ==================================================================
C ==================================================================
SUBROUTINE vecrotx(angle,vec1,vec2)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION vec1(3), vec2(3)
REAL*8 R(3,3)
C ==--------------------------------------------------------------==
DO i=1,3
DO j=1,3
R(i,j)=0.0D0
END DO
R(i,i)=1.0D0
vec2(i)=0.0D0
END DO
R(2,2)=+cos(angle)
R(3,3)=+cos(angle)
R(2,3)=-sin(angle)
R(3,2)=+sin(angle)
DO i=1,3
DO j=1,3
vec2(I)=vec2(I)+R(I,J)*vec1(J)
END DO
END DO
C ==--------------------------------------------------------------==
RETURN
END
C ==================================================================
C
C ==================================================================
SUBROUTINE vecrotvec(vec1,angle,vec2,vec3)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION vec1(3), vec2(3), vec3(3)
REAL*8 dummy(3,6)
REAL*8 a1,a2,a3,pi
C ==--------------------------------------------------------------==
pi=acos(-1.0D0)
a1=(pi/2.0D0)-atan2(vec1(2),vec1(1)) ! first y than x
call vecrotz(a1,vec1,dummy(1,1)) ! rotate in the yz plane
call vecrotz(a1,vec2,dummy(1,2))
a2=(pi/2.0D0)-atan2(dummy(3,1),dummy(2,1))
call vecrotx(a2,dummy(1,1),dummy(1,3)) ! put on the z axis
call vecrotx(a2,dummy(1,2),dummy(1,4))
! rotate now around z for the angle
call vecrotz(angle,dummy(1,4),dummy(1,5))
! and go back
call vecrotx(-a2,dummy(1,5),dummy(1,6))
call vecrotz(-a1,dummy(1,6),vec3)
C ==--------------------------------------------------------------==
RETURN
END
C ==================================================================
More information about the CPMD-list
mailing list