1      subroutine planky(npts1,npts2,npts3,keywrd,defau)
2      implicit double precision (a-h,o-z)
3      parameter (numatm=2000)
4      parameter (mxvalc=10)
5      logical defau
6      logical keyi,keyr,keyiv,keyrv,keyirv
7      character*(*) keywrd
8      common /rdwr/   iun1,iun2,iun3,iun4,iun5
9      common /plane/  px, py, pz, cx, cy, cz, r(3),v1(3),v2(3),iplat
10      common /coord / xyz(3,numatm)
11      common /moldat/ natoms, norbs, nelecs, nat(numatm)
12      logical oscal,ostep,fine
13      common /plspec/ step,stepf,scale,scalf,cut,pmin,pmax,one,
14     &                oscal,ostep,fine
15      common /srfhlp/edge,ctval(mxvalc),pxyz(3),nvalc,nspts,istyp
16      common /grdhlp/ mx3d,mx3d2
17      dimension d1(3),d2(3),c12(3)
18
19      todeg = 45.0d0 / datan(1.0d0)
20      iflag = 0
21      imxflg = 0
22
23      if (.not.keyrv(keywrd,'CENTER',px,py,pz)) then
24          if (keyi(keywrd,'CENTER',i)) then
25             if (i.gt.natoms) then
26                call inferr('CENTER: Atom Nr. > Nr. of Atoms !',1)
27                return
28             endif
29             iplat = 0
30             iflag = 1
31             px = xyz(1,i)
32             py = xyz(2,i)
33             pz = xyz(3,i)
34          endif
35      else
36          iplat = 0
37          iflag = 1
38      endif
39
40      if (.not.keyrv(keywrd,'LINE',cx,cy,cz)) then
41          if (keyi(keywrd,'LINE',i)) then
42             if (i.gt.natoms) then
43                call inferr('LINE: Atom Nr. > Nr. of Atoms !',1)
44                return
45             endif
46             iplat = 0
47             iflag = 1
48             cx = xyz(1,i) - px
49             cy = xyz(2,i) - py
50             cz = xyz(3,i) - pz
51          endif
52      else
53          iplat = 0
54          iflag = 1
55      endif
56
57      if (keyiv(keywrd,'PLANE',iat1,iat2,iat3)) then
58          if (iflag.eq.1) then
59             call inferr('CENTER/LINE and Plane EXCLUSIVE!',1)
60             return
61          endif
62          if((iat1.gt.natoms).or.(iat2.gt.natoms).or.(iat3.gt.natoms))
63     &    then
64             call inferr('PLANE: Atom nr. > number of atoms !',1)
65             return
66          endif
67          call parpla(iat1,iat2,iat3,istat)
68          if (istat.eq.1) print*,'Three atoms lie on a line !'
69          iflag = 2
70      endif
71
72      if (iflag.eq.0.and.defau) call defpc
73
74      if (keyirv(keywrd,'ROT',iata,iatb,degree)) then
75          if (iflag.eq.1.or.iflag.eq.0) then
76            if (iflag.eq.1)
77     &         call inferr('Center/Line and Rot are EXCLUSIVE',1)
78            if (iflag.eq.0)
79     &         call inferr('ROT used before PLANE !',1)
80            return
81          endif
82          if (((iata.ne.iat1).and.(iata.ne.iat2).and.(iata.ne.iat3))
83     &    .or.((iatb.ne.iat1).and.(iatb.ne.iat2).and.(iatb.ne.iat3)))
84     &    then
85             call inferr('ROT: atom used not used with Plane !',1)
86             return
87          endif
88          radian = degree / todeg
89          do i=1,3
90            d1(i) = (xyz(i,iatb) - xyz(i,iata))
91          end do
92          d2(1) = cx
93          d2(2) = cy
94          d2(3) = cz
95          px = xyz(1,iata) + 0.5d0*d1(1)
96          py = xyz(2,iata) + 0.5d0*d1(2)
97          px = xyz(3,iata) + 0.5d0*d1(3)
98          d2norm = vlen(d2)
99          do i=1,3
100            d2(i) = d2(i)/d2norm
101          end do
102          call crprod(d1,d2,c12)
103          c12nrm = vlen(c12)
104          do i=1,3
105            c12(i) = c12(i)/c12nrm
106          end do
107          cx = dcos(radian)*d2(1) - dsin(radian)*c12(1)
108          cy = dcos(radian)*d2(2) - dsin(radian)*c12(2)
109          cz = dcos(radian)*d2(3) - dsin(radian)*c12(3)
110          iplat = 0
111      endif
112
113      if (keyr(keywrd,'LIFT',rlift)) call dolift(rlift)
114
115      if (.not.keyr(keywrd,'EDGE',r(1))) then
116         if (defau) call defrad(.false.)
117      else
118         r(2) = r(1)
119         r(3) = r(1)
120      endif
121
122      if (.not.keyr(keywrd,'EDZ',r(3))) then
123         r(3) = r(1)
124      endif
125
126      if (.not.keyr(keywrd,'EDY',r(2))) then
127         r(2) = r(1)
128      endif
129
130      if (.not.keyr(keywrd,'EDX',r(1))) idum = 1
131
132      if (r(1).lt.0.01d0) then
133         r(1) = 3.0
134         r(2) = r(1)
135         r(3) = r(1)
136      endif
137
138      if (keyi(keywrd,'NPTSX',npts1)) then
139         if (npts1.gt.mx3d) imxflg = 1
140      endif
141      if (keyi(keywrd,'NPTSY',npts2)) then
142         if (npts2.gt.mx3d) imxflg = 1
143      endif
144      if (keyi(keywrd,'NPTSZ',npts3)) then
145         if (npts3.gt.mx3d) imxflg = 1
146      endif
147
148      if (imxflg.eq.1) then
149         nptsmx = npts1
150         if (npts2.gt.nptsmx) nptsmx = npts2
151         if (npts3.gt.nptsmx) nptsmx = npts3
152         call allgrd(nptsmx)
153         if (nptsmx.gt.mx3d) then
154             call inferr('Exceeding maximum points !',1)
155             if (npts1.gt.mx3d) npts1 = mx3d
156             if (npts2.gt.mx3d) npts2 = mx3d
157             if (npts3.gt.mx3d) npts3 = mx3d
158         endif
159      endif
160
161      oneold = one
162      if (.not.keyr(keywrd,'PHASE',one)) then
163         if (index(keywrd,'PHASE').ne.0) one = -1.d0*one
164      else
165         if (one.eq.0.0d0) one = -1.d0*oneold
166      endif
167
168      if (index(keywrd,'ALIGN').ne.0) then
169         px = pxyz(1)
170         py = pxyz(2)
171         pz = pxyz(3)
172         iplat = 0
173      endif
174
175      write(iun3,10) px, py, pz, cx, cy, cz, r(1), r(2), r(3)
17610    format(9x,'  CENTER OF GRAPH           AXIS OF GRAPH       RADII',
177     1/9X,      'PX      PY      PZ        CX      CY      CZ',
178     24X,      'EDX   EDY   EDY',
179     3/,5X,3F8.3,F10.3,2F8.3,3F6.1,///)
180
181      return
182      end
183
184      logical function keyi(keyl,keyw,ival)
185      implicit double precision (a-h,o-z)
186      character*(*) keyl,keyw
187
188      keyi = .false.
189      l = linlen(keyw)
190      i = index(keyl,keyw)
191      if (i.ne.0) then
192          if ((i.eq.1.or.keyl(i-1:i-1).eq.' ').and.
193     &         keyl(i+l:i+l).eq.'=') then
194             ival = reada(keyl,i,len(keyl))
195             keyi = .true.
196          endif
197      endif
198
199      return
200      end
201
202      logical function keyr(keyl,keyw,rval)
203      implicit double precision (a-h,o-z)
204      character*(*) keyl,keyw
205
206      keyr = .false.
207      i = index(keyl,keyw)
208      l = linlen(keyw)
209      if (i.ne.0) then
210          if ((i.eq.1.or.keyl(i-1:i-1).eq.' ').and.
211     &         keyl(i+l:i+l).eq.'=') then
212             rval = reada(keyl,i,len(keyl))
213             keyr = .true.
214          endif
215      endif
216
217      return
218      end
219
220      logical function keyrv(keyl,keyw,r1,r2,r3)
221      implicit double precision (a-h,o-z)
222      character*(*) keyl,keyw
223
224      keyrv = .false.
225      kend = len(keyl)
226
227      i = index(keyl,keyw)
228      if (i.eq.0) return
229      i = i + len(keyw)
230      j = index(keyl(i:kend),'(')
231      if (j.eq.0) return
232      do k=0,j-2
233         if (keyl(i+k:i+k).ne.' '.and.keyl(i+k:i+k).ne.'=') return
234      end do
235      i = i + j
236      r1 = reada(keyl,i,len(keyl))
237      j = index(keyl(i:kend),',')
238      if (j.eq.0) return
239      i = i + j
240      r2 = reada(keyl,i,len(keyl))
241      j = index(keyl(i:kend),',')
242      if (j.eq.0) return
243      i = i + j
244      r3 = reada(keyl,i,len(keyl))
245
246      keyrv = .true.
247      return
248      end
249
250      integer function keyr3v(keyl,keyw,r,n)
251      implicit double precision (a-h,o-z)
252      character*(*) keyl,keyw
253      dimension r(*)
254
255      keyr3v = 0
256      kend = len(keyl)
257
258      i = index(keyl,keyw)
259      if (i.eq.0) return
260      i = i + len(keyw)
261      j = index(keyl(i:kend),'(')
262      if (j.eq.0) return
263      do k=0,j-2
264         if (keyl(i+k:i+k).ne.' '.and.keyl(i+k:i+k).ne.'=') return
265      end do
266
267      do l=1,n
268         i = i + j
269         r(l) = reada(keyl,i,len(keyl))
270         keyr3v = keyr3v + 1
271         j = index(keyl(i:kend),',')
272         if (j.eq.0) return
273      end do
274
275      return
276      end
277
278      logical function keyiv(keyl,keyw,i1,i2,i3)
279      implicit double precision (a-h,o-z)
280      character*(*) keyl,keyw
281
282      keyiv = .false.
283      kend = len(keyl)
284
285      i = index(keyl,keyw)
286      if (i.eq.0) return
287      i = i + len(keyw)
288      j = index(keyl(i:kend),'(')
289      if (j.eq.0) return
290      do k=0,j-2
291         if (keyl(i+k:i+k).ne.' '.and.keyl(i+k:i+k).ne.'=') return
292      end do
293      i = i + j
294      i1 = reada(keyl,i,len(keyl))
295      j = index(keyl(i:kend),',')
296      if (j.eq.0) return
297      i = i + j
298      i2 = reada(keyl,i,len(keyl))
299      j = index(keyl(i:kend),',')
300      if (j.eq.0) return
301      i = i + j
302      i3 = reada(keyl,i,len(keyl))
303
304      keyiv = .true.
305      return
306      end
307
308      logical function keyirv(keyl,keyw,i1,i2,r)
309      implicit double precision (a-h,o-z)
310      character*(*) keyl,keyw
311
312      keyirv = .false.
313      kend = len(keyl)
314
315      i = index(keyl,keyw)
316      if (i.eq.0) return
317      i = i + len(keyw)
318      j = index(keyl(i:kend),'(')
319      if (j.eq.0) return
320      do k=0,j-2
321         if (keyl(i+k:i+k).ne.' '.and.keyl(i+k:i+k).ne.'=') return
322      end do
323      i = i + j
324      i1 = reada(keyl,i,len(keyl))
325      j = index(keyl(i:kend),',')
326      if (j.eq.0) return
327      i = i + j
328      i2 = reada(keyl,i,len(keyl))
329      j = index(keyl(i:kend),',')
330      if (j.eq.0) return
331      i = i + j
332      r = reada(keyl,i,len(keyl))
333
334      keyirv = .true.
335      return
336      end
337