1c
2c $Id$
3c
4      subroutine nucdd_calc(c,zan,nat,dde)
5c
6c  This is a modified HONDO routine (hnd_nucdx) to calculate
7c  the second derivative contribution of the nuclear term to
8c  the hessian (dde).
9c
10      implicit none
11c
12#include "global.fh"
13c
14      integer nat                       ! [input] number of atoms
15      double precision zan(nat)         ! [input] nuclear charges
16      double precision c(3,nat)         ! [input] cartesian coordinates
17      double precision dde(3,nat,3,nat) ! [input/output] hessian matrix
18c
19      integer iat, jat, ixyz, myid, nprocs
20      double precision rij, dum, dum1, zero, three
21c
22      data zero   /0.0d+00/
23      data three  /3.0d+00/
24c
25      if(nat.eq.1) go to 200
26c
27c     ----- calculate forces -----
28c
29      myid = ga_nodeid()  ! set up static parallel constructs
30      nprocs = ga_nnodes()
31c
32      do 140 iat=myid+2,nat,nprocs   ! go parallel
33      do 140 jat=1,iat-1
34c
35        rij=zero
36        do 110 ixyz=1,3
37          rij=rij+(c(ixyz,iat)-c(ixyz,jat))**2
38  110   continue
39        dum = zan(iat)*zan(jat)/(rij*sqrt(rij))
40c
41c     ----- ddedii , ddedij , ddedjj -----
42c
43        dum=dum/rij
44        do 130 ixyz=1,3
45          dum1=(three*(c(ixyz,iat)-c(ixyz,jat))**2-rij)*dum
46          dde(ixyz,iat,ixyz,iat)=dde(ixyz,iat,ixyz,iat)+dum1
47          dde(ixyz,iat,ixyz,jat)=dde(ixyz,iat,ixyz,jat)-dum1
48          dde(ixyz,jat,ixyz,iat)=dde(ixyz,jat,ixyz,iat)-dum1
49          dde(ixyz,jat,ixyz,jat)=dde(ixyz,jat,ixyz,jat)+dum1
50  130   continue
51        dum1=(three*(c(1,iat)-c(1,jat))*(c(2,iat)-c(2,jat)))*dum
52        dde(1,iat,2,iat)=dde(1,iat,2,iat)+dum1
53        dde(2,iat,1,iat)=dde(2,iat,1,iat)+dum1
54        dde(1,iat,2,jat)=dde(1,iat,2,jat)-dum1
55        dde(2,iat,1,jat)=dde(2,iat,1,jat)-dum1
56        dde(1,jat,2,jat)=dde(1,jat,2,jat)+dum1
57        dde(2,jat,1,jat)=dde(2,jat,1,jat)+dum1
58c
59        if (iat.ne.jat) then
60          dde(1,jat,2,iat)=dde(1,jat,2,iat)-dum1
61          dde(2,jat,1,iat)=dde(2,jat,1,iat)-dum1
62        endif
63c
64        dum1=(three*(c(1,iat)-c(1,jat))*(c(3,iat)-c(3,jat)))*dum
65        dde(1,iat,3,iat)=dde(1,iat,3,iat)+dum1
66        dde(3,iat,1,iat)=dde(3,iat,1,iat)+dum1
67        dde(1,iat,3,jat)=dde(1,iat,3,jat)-dum1
68        dde(3,iat,1,jat)=dde(3,iat,1,jat)-dum1
69        dde(1,jat,3,jat)=dde(1,jat,3,jat)+dum1
70        dde(3,jat,1,jat)=dde(3,jat,1,jat)+dum1
71c
72        if (iat.ne.jat) then
73          dde(1,jat,3,iat)=dde(1,jat,3,iat)-dum1
74          dde(3,jat,1,iat)=dde(3,jat,1,iat)-dum1
75        endif
76c
77        dum1=(three*(c(2,iat)-c(2,jat))*(c(3,iat)-c(3,jat)))*dum
78        dde(2,iat,3,iat)=dde(2,iat,3,iat)+dum1
79        dde(3,iat,2,iat)=dde(3,iat,2,iat)+dum1
80        dde(2,iat,3,jat)=dde(2,iat,3,jat)-dum1
81        dde(3,iat,2,jat)=dde(3,iat,2,jat)-dum1
82        dde(2,jat,3,jat)=dde(2,jat,3,jat)+dum1
83        dde(3,jat,2,jat)=dde(3,jat,2,jat)+dum1
84c
85        if (iat.ne.jat) then
86          dde(2,jat,3,iat)=dde(2,jat,3,iat)-dum1
87          dde(3,jat,2,iat)=dde(3,jat,2,iat)-dum1
88        endif
89c
90  140 continue
91c
92  200 continue
93c
94      return
95      end
96