forked from hornos/lev00
-
Notifications
You must be signed in to change notification settings - Fork 0
/
device.f90
707 lines (680 loc) · 24.6 KB
/
device.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
subroutine check_atoms(nat,good)
use atoms
implicit none
!.............................................................
! checks if atom nat is equivalent to any of k=1,...,nat-1
!.............................................................
real*8 r(3)
integer nat,k
logical good,ask_equiv
if(nat.gt.1) then
do k=1,nat-1
r(1)=TI(1,k)-TI(1,nat)
r(2)=TI(2,k)-TI(2,nat)
r(3)=TI(3,k)-TI(3,nat)
if(ask_equiv(r,BCELL,3,tiny_equiv)) then
write(*,'(a17,i2,a27,i2)')'ERROR! Your atom ', &
nat,' is equivalent to atom ',k
good=.false.
end if
end do
end if
return
end subroutine check_atoms
subroutine ask(iX,iY,iZ,Point,Radius,DIRC,iask)
!....................................................................
! Replies with:
! iask='y' - if the grid point (iX,iY,iZ) is inside a sphere
! of the Radius drawn around the Point;
! iask='n' - if not.
!....................................................................
! DIRC - direct lattice vectors
!....................................................................
use param
implicit none
character iask*1
integer i,j,k,iX,iY,iZ
real*8 Point(3),r(3),x(3),DIRC(3,3),Radius,dist
!
! _______ get the point (iX,iY,iZ) in Cartesian coordinates as folows:
r(1) = iX*DIRC(1,1)/NGX + iY*DIRC(2,1)/NGY + iZ*DIRC(3,1)/NGZ
r(2) = iX*DIRC(1,2)/NGX + iY*DIRC(2,2)/NGY + iZ*DIRC(3,2)/NGZ
r(3) = iX*DIRC(1,3)/NGX + iY*DIRC(2,3)/NGY + iZ*DIRC(3,3)/NGZ
!
!_______ a distance, dist, between r() and 27 equivalent spheres
! with the centers put in x(), is calculated. The latter are built
! from Point (the center of the input sphere) by adding the smallest
! 27 direct lattice vectors:
!
do i=-1,1
do j=-1,1
do k=-1,1
x(1) = i*DIRC(1,1)+j*DIRC(2,1)+k*DIRC(3,1) + Point(1)
x(2) = i*DIRC(1,2)+j*DIRC(2,2)+k*DIRC(3,2) + Point(2)
x(3) = i*DIRC(1,3)+j*DIRC(2,3)+k*DIRC(3,3) + Point(3)
dist=sqrt((r(1)-x(1))**2+(r(2)-x(2))**2+(r(3)-x(3))**2)
if(dist.lt.Radius) then
iask='y'
return
end if
end do
end do
end do
iask='n'
return
end subroutine ask
subroutine ask_box(iX,iY,iZ,Center,SIDES,RECIP,DIRC,iask)
!.........................................................................
! Replies with:
! iask='y' - if the grid point (iX,iY,iZ) is inside a parallelepiped
! (the box) drawn around the Center with sides along the
! directions of the vectors SIDES
! iask='n' - if not.
!.........................................................................
! DIRC(#,xyz) - direct lattice vectors
! SIDES(#,xyz) - directions of the box sides
! RECIP(#,xyz) - "reciprocal" vectors associated with SIDES, i.e.
! RECIP(#1)*SIDES(#2)=delta(#1,#2) (dot product)
!.........................................................................
use param
implicit none
character iask*1
integer iX,iY,iZ
real*8 Center(3),DIRC(3,3),r(3),SIDES(3,3),RECIP(3,3)
!
!_______ get the point (iX,iY,iZ) in Cartesian coordinates as folows:
r(1) = iX*DIRC(1,1)/NGX + iY*DIRC(2,1)/NGY + iZ*DIRC(3,1)/NGZ
r(2) = iX*DIRC(1,2)/NGX + iY*DIRC(2,2)/NGY + iZ*DIRC(3,2)/NGZ
r(3) = iX*DIRC(1,3)/NGX + iY*DIRC(2,3)/NGY + iZ*DIRC(3,3)/NGZ
!
!_______ ask if the point r() is inside the box
call ask_box1(r,Center,SIDES,RECIP,DIRC,iask)
end subroutine ask_box
subroutine ask_box1(r,Center,SIDES,RECIP,DIRC,iask)
!.........................................................................
! Replies with:
! iask='y' - if the Cartesian grid point r() is inside a parallelepiped
! (the box) drawn around the Center with sides along the
! directions of the vectors SIDES
! iask='n' - if not.
!.........................................................................
! DIRC(#,xyz) - direct lattice vectors
! SIDES(#,xyz) - directions of the box sides
! RECIP(#,xyz) - "reciprocal" vectors associated with SIDES, i.e.
! RECIP(#1)*SIDES(#2)=delta(#1,#2) (dot product)
!.........................................................................
use param
implicit none
character iask*1
real*8,parameter :: tiny = 0.0001
real*8 Center(3),DIRC(3,3),r(3),x(3),corner(3),SIDES(3,3)
real*8 rins(3),RECIP(3,3),d0,d1,c1,c2,c3
integer i,j,k
!
!_______ the vector to the 1st corner of the box
!
corner(1)=Center(1)-0.5*(SIDES(1,1)+SIDES(2,1)+SIDES(3,1))
corner(2)=Center(2)-0.5*(SIDES(1,2)+SIDES(2,2)+SIDES(3,2))
corner(3)=Center(3)-0.5*(SIDES(1,3)+SIDES(2,3)+SIDES(3,3))
!
!_______ part of r() inside the box
!
rins(1)=r(1)-corner(1)
rins(2)=r(2)-corner(2)
rins(3)=r(3)-corner(3)
!
!_______ The vector rins() is inside the box fixed by vectors SIDES if
! in the expansion of rins() wrt SIDES all 3 coefficients are
! between 0.0 and 1.0
!_______ We run over all possible 27 images of the box to check if the
! point rins() gets inside at least one of those.
!
d0=-tiny
d1=1.0+tiny
do i=-1,1
do j=-1,1
do k=-1,1
x(1) = i*DIRC(1,1)+j*DIRC(2,1)+k*DIRC(3,1) + rins(1)
x(2) = i*DIRC(1,2)+j*DIRC(2,2)+k*DIRC(3,2) + rins(2)
x(3) = i*DIRC(1,3)+j*DIRC(2,3)+k*DIRC(3,3) + rins(3)
c1=x(1)*RECIP(1,1)+x(2)*RECIP(1,2)+x(3)*RECIP(1,3)
if(c1.ge.d0 .and. c1.le.d1) then
c2=x(1)*RECIP(2,1)+x(2)*RECIP(2,2)+x(3)*RECIP(2,3)
if(c2.ge.d0 .and. c2.le.d1) then
c3=x(1)*RECIP(3,1)+x(2)*RECIP(3,2)+x(3)*RECIP(3,3)
if(c3.ge.d0 .and. c3.le.d1) then
iask='y'
return
end if
end if
end if
end do
end do
end do
iask='n'
end subroutine ask_box1
subroutine ask_box2(r,Center,SIDES,RECIP,iask)
!.........................................................................
! Replies with:
! iask='y' - if the grid point r() is inside a parallelepiped
! (the box) drawn around the Center with sides along the
! directions of the vectors SIDES
! iask='n' - if not.
!.........................................................................
! SIDES(#,xyz) - directions of the box sides
! RECIP(#,xyz) - "reciprocal" vectors associated with SIDES, i.e.
! RECIP(#1)*SIDES(#2)=delta(#1,#2) (dot product)
!.........................................................................
! the only difference with ask_box1() is that images are not used.
!.........................................................................
use param
implicit none
character iask*1
real*8,parameter :: tiny=0.0001
real*8 Center(3),r(3),x(3),corner(3),SIDES(3,3),RECIP(3,3),d0,d1,c1,c2,c3
!
!_______ the vector to the 1st corner of the box
!
corner(1)=Center(1)-0.5*(SIDES(1,1)+SIDES(2,1)+SIDES(3,1))
corner(2)=Center(2)-0.5*(SIDES(1,2)+SIDES(2,2)+SIDES(3,2))
corner(3)=Center(3)-0.5*(SIDES(1,3)+SIDES(2,3)+SIDES(3,3))
!
!_______ part of r() inside the box
!
x(1)=r(1)-corner(1)
x(2)=r(2)-corner(2)
x(3)=r(3)-corner(3)
!
!_______ The vector x() is inside the box fixed by vectors SIDES if
! in the expansion of x() wrt SIDES all 3 coefficients are
! between 0.0 and 1.0
d0=-tiny
d1=1.0+tiny
c1=x(1)*RECIP(1,1)+x(2)*RECIP(1,2)+x(3)*RECIP(1,3)
if(c1.ge.d0 .and. c1.le.d1) then
c2=x(1)*RECIP(2,1)+x(2)*RECIP(2,2)+x(3)*RECIP(2,3)
if(c2.ge.d0 .and. c2.le.d1) then
c3=x(1)*RECIP(3,1)+x(2)*RECIP(3,2)+x(3)*RECIP(3,3)
if(c3.ge.d0 .and. c3.le.d1) then
iask='y'
return
end if
end if
end if
iask='n'
end subroutine ask_box2
subroutine show_atoms()
use param
use atoms
implicit none
integer k,i,j,j1,i1
real*8 x(3)
k=0
do i=1,NSPEC
write(*,'(3a,i5)')'.... Species <',Species(i),'> ......> ',i
write(*,'(a)')' Tot # |---------- Fractional ---------|' &
//'------------- Cartesian ---------|'
do j=1,NspN(i)
k=k+1
do i1=1,3
x(i1)=0.0
do j1=1,3
x(i1)=x(i1)+TI(j1,k)*BCELL(i1,j1)
end do
end do
write(*,'(2(x,i3),x,3(x,f10.5),2x,3(x,f10.5),3x,a,4x,a)') &
k,j,(x(i1),i1=1,3),(TI(j1,k),j1=1,3)
end do
end do
write(*,*)'Hit ENTER when done ...'
read(*,*)
end subroutine show_atoms
subroutine show_atoms_tagged(listat,ngroup)
use param
use atoms
implicit none
logical listat(NIONS,ngroup)
character cha
integer k,i,j,j1,ngroup
real*8 a1,a2,a3
k=0
do i=1,NSPEC
write(*,'(3a,i5)')'.... Species <',Species(i),'> ......> ',i
write(*,'(a)')' Tot # |---------- Fractional ---------|' &
//'------------- Cartesian ---------| Tag'
do j=1,NspN(i)
k=k+1
a1=TI(1,k)*BCELL(1,1)+TI(2,k)*BCELL(1,2)+TI(3,k)*BCELL(1,3)
a2=TI(1,k)*BCELL(2,1)+TI(2,k)*BCELL(2,2)+TI(3,k)*BCELL(2,3)
a3=TI(1,k)*BCELL(3,1)+TI(2,k)*BCELL(3,2)+TI(3,k)*BCELL(3,3)
cha=' '
do j1=1,ngroup
if(listat(k,j1)) cha='Y'
end do
write(*,'(2(x,i3),x,3(x,f10.5),2x,3(x,f10.5),3x,a)') &
k,j,a1,a2,a3,(TI(j1,k),j1=1,3),cha
end do
end do
write(*,*)'Hit ENTER when done ...'
read(*,*)
end subroutine show_atoms_tagged
subroutine show_atoms_in_plane()
!...............................................................................
! Atomic positions with respect to the plane specified by the normal
! (vers1x,vers1y,vers1z) and the cnetral point (aCENTX,aCENTY,aCENTZ)
! are given for all atoms using the local coordinate system with:
! X axis = (vers2x,vers2y,vers2z)
! Y axis = (vers3x,vers3y,vers3z)
!...............................................................................
use param
use atoms
use menu
implicit none
real*8,parameter :: tiny=0.0001
real*8 x(3),cnt(3)
integer k,i,j,i1,j1
character cha*10
!
cnt(1)=aCENTX
cnt(2)=aCENTY
cnt(3)=aCENTZ
!
k=0
do i=1,NSPEC
write(*,'(3a,i5)')'.... Species <',Species(i),'> ......> ',i
write(*,'(a)')' Tot # |---------- Cartesian ----------|' &
//' X-plane Y-plane out-of-plane'
do j=1,NspN(i)
k=k+1
!______________ position x() of atom k in the coordinate system fixed to the
! plane central point (aCENTX,aCENTY,aCENTZ) and axes X,Y,and Z=normal
! (x(3) shows out-of-plane distance)
!
x(1)=(TI(1,k)-aCENTX)*vers2x+(TI(2,k)-aCENTY)*vers2y+ &
(TI(3,k)-aCENTZ)*vers2z
x(2)=(TI(1,k)-aCENTX)*vers3x+(TI(2,k)-aCENTY)*vers3y+ &
(TI(3,k)-aCENTZ)*vers3z
x(3)=(TI(1,k)-aCENTX)*vers1x+(TI(2,k)-aCENTY)*vers1y+ &
(TI(3,k)-aCENTZ)*vers1z
if(abs(x(3)).lt.tiny) then
cha=' NO '
else
write(cha,'(f10.5)') x(3)
end if
write(*,'(2(x,i3),x,3(x,f10.5),2x,2(x,f10.5),x,a)') &
k,j,(TI(j1,k),j1=1,3),(x(i1),i1=1,2),cha
end do
end do
write(*,*)'Hit ENTER when done ...'
read(*,*)
end subroutine show_atoms_in_plane
subroutine at_in_box(Cent,Sides,Rec,Pos,iType,N_at,nt,Err)
!.................................................................
! Positions Pos() and species types iType() of all atoms fitted
! into the box (with Cent, Sides, Rec) are calculated.
!.................................................................
! nt - number of atoms found in the box
!.................................................................
use param
use atoms
implicit none
real*8 Pos(3,N_at),Cent(3),Sides(3,3),Rec(3,3),x(3)
integer iType(N_at),i,j,k,i1,j1,nat,nt,iErr,N_at
character iask
logical Err
Err=.false.
nat=0
nt=0
DO i1=1,NSPEC
DO j1=1,NspN(i1)
nat=nat+1
do i=-1,1
do j=-1,1
do k=-1,1
x(1)=i*DIRC(1,1)+j*DIRC(2,1)+k*DIRC(3,1)+TI(1,nat)
x(2)=i*DIRC(1,2)+j*DIRC(2,2)+k*DIRC(3,2)+TI(2,nat)
x(3)=i*DIRC(1,3)+j*DIRC(2,3)+k*DIRC(3,3)+TI(3,nat)
call ask_box2(x,Cent,Sides,Rec,iask)
if(iask.eq.'y') then
nt=nt+1
if(nt.eq.N_at) then
write(*,*) &
'ERROR! # of atoms in the box > ',N_at
Err=.true.
return
end if
Pos(1,nt)=x(1)
Pos(2,nt)=x(2)
Pos(3,nt)=x(3)
iType(nt)=i1
end if
end do
end do
end do
END DO
END DO
write(*,*)'Number of atoms in the box found = ',nt
end subroutine at_in_box
subroutine reducn(x,DIRC,BCELL)
!..................................................................
! It returns the point x() to the 0-th unit cell by extracting
! integer parts from the coefficients in the decomposition of x()
! with respect to DIRC and by making them positive
!..................................................................
implicit none
real*8 x(3),DIRC(3,3),BCELL(3,3),r(3)
real*8,parameter :: tiny=0.000001
integer i
!
!________ get coefficients r() in the decomposition of x() with
! respect to DIRC;
! These coefficients are iX/NGX, iY/NGY, iZ/NGZ, where (iX,iY,iZ)
! might be larger than NGX, NGY,NGZ and noninteger; this depends on
! the point x() which has nothing to do with the grid.
!
r(1)=BCELL(1,1)*x(1)+BCELL(1,2)*x(2)+BCELL(1,3)*x(3)
r(2)=BCELL(2,1)*x(1)+BCELL(2,2)*x(2)+BCELL(2,3)*x(3)
r(3)=BCELL(3,1)*x(1)+BCELL(3,2)*x(2)+BCELL(3,3)*x(3)
!
!________ extract the integer part from the coefficients and make
! the coefficients positive if they are not
!
do i=1,3
if(abs(r(i)-nint(r(i))).lt.tiny) r(i)=nint(r(i))
r(i)=r(i)-int(r(i))+1
r(i)=r(i)-int(r(i))
end do
!
!________ return the point r() to the Cartesian representation, x()
!
x(1)=r(1)*DIRC(1,1)+r(2)*DIRC(2,1)+r(3)*DIRC(3,1)
x(2)=r(1)*DIRC(1,2)+r(2)*DIRC(2,2)+r(3)*DIRC(3,2)
x(3)=r(1)*DIRC(1,3)+r(2)*DIRC(2,3)+r(3)*DIRC(3,3)
end subroutine reducn
subroutine normalize(x,y,z)
implicit none
real*8 znorm,x,y,z
znorm=sqrt(x*x+y*y+z*z)
x=x/znorm
y=y/znorm
z=z/znorm
end subroutine normalize
subroutine interpolate(R,BCELL,denval,grid)
!.......................................................
! This subroutine does the internal interpolation
! of the point R within the eight cube vertices of the grid.
!.......................................................
! Every point in the grid R(i,j,k) is defined as:
!
! R(i,j,k)=DIRC(1)*i/NGX+DIRC(2)*j/NGY+DIRC(3)*k/NGZ
!
! where i = 0,...,NGX-1; j = 0,...,NGY-1;
! k = 0,...,NGZ-1 - are used to numerate
! the grid.
!.......................................................
use param
implicit none
real*8 R(3), x(3), BCELL(3,3)
integer intx(3), cube(2,2,2,3),i1,i2,i3
real*8 grid(NGX,NGY,NGZ), frac(3),t5,t6,denval,t1,t2,t3,t4
integer,dimension(3) :: nmesh
nmesh(1)=NGX ; nmesh(2)=NGY ; nmesh(3)=NGZ
!
!________ get coefficients x() in the decomposition of r() with
! respect to DIRC;
! These coefficients are iX/NGX, iY/NGY, iZ/NGZ, where (iX,iY,iZ)
! are exactly between 0 and NGX-1, 0 and NGY-1, etc. because of
! reducn() called previously
!
x(1)=BCELL(1,1)*r(1)+BCELL(1,2)*r(2)+BCELL(1,3)*r(3)
x(2)=BCELL(2,1)*r(1)+BCELL(2,2)*r(2)+BCELL(2,3)*r(3)
x(3)=BCELL(3,1)*r(1)+BCELL(3,2)*r(2)+BCELL(3,3)*r(3)
!
!....... define position of x() in the net of grid points
!
do i1=1,3
x(i1)=nmesh(i1)*x(i1)
frac(i1)=x(i1)-int(x(i1))
intx(i1)=int(x(i1))
end do
!
!___________ put x() inside a cube with vertices on the net and get
! its numbers in cube(i1,i2,i3,component), where (i1,i2,i3) give all
! possible 8 vertices of the cube (because of grid(iX+1,iY+1,iZ+1)
! defined for indices 1,...NGX, etc., the array cube gets +1):
!
do i1=1,2
do i2=1,2
do i3=1,2
cube(i1,i2,i3,1)=intx(1)+i1
cube(i1,i2,i3,2)=intx(2)+i2
cube(i1,i2,i3,3)=intx(3)+i3
if(cube(i1,i2,i3,1).eq.(NGX+1)) cube(i1,i2,i3,1)=1
if(cube(i1,i2,i3,2).eq.(NGY+1)) cube(i1,i2,i3,2)=1
if(cube(i1,i2,i3,3).eq.(NGZ+1)) cube(i1,i2,i3,3)=1
end do
end do
end do
!
!____________ interpolate the density in x() from that on the cube
! vectices
!
call lineint(cube,grid,1,1,1,2,1,1,frac(1),t1)
call lineint(cube,grid,1,2,1,2,2,1,frac(1),t2)
call lineint(cube,grid,1,1,2,2,1,2,frac(1),t3)
call lineint(cube,grid,1,2,2,2,2,2,frac(1),t4)
t5=t1+frac(2)*(t2-t1)
t6=t3+frac(2)*(t4-t3)
denval=t5+frac(3)*(t6-t5)
end subroutine interpolate
subroutine lineint(cube,grid,p1,p2,p3,q1,q2,q3,frac,val)
!.......................................................
! A small subroutine to do linear interpolation between
! two cube vertices
!.......................................................
use param
implicit none
integer cube(2,2,2,3),p1,p2,p3,q1,q2,q3
real*8 grid(NGX,NGY,NGZ),val,t1,t2,frac
t1=grid(cube(p1,p2,p3,1),cube(p1,p2,p3,2),cube(p1,p2,p3,3))
t2=grid(cube(q1,q2,q3,1),cube(q1,q2,q3,2),cube(q1,q2,q3,3))
val=t1+frac*(t2-t1)
end subroutine lineint
subroutine vector3(DIRC)
!.....................................................................
! The plane case:
! Choses a normal vector (vers1x,vers1y,vers1z) to the plane as well
! as two orthonormal vectors (vers2x,...) and (vers3x,...) in the
! plane.
!.....................................................................
! 'icase' gives the method used to specify the plane:
! icase = 1 - by a normal vector;
! 2 - by 3 points lying in the plane.
! angstr - gives a method of specifying coordinates
!.....................................................................
use menu
implicit none
real*8 Rbb(3),Rcc(3),DIRC(3,3),a
character answer
real*8,parameter :: tiny=0.00001
!
!____________ plane: a vector along the normal to the plane
!
!_________________ if the plane is chosen by the normal vector
if(icase.eq.1) then
7 write(*,*)'Give a vector (x,y,z) orthogonal to your plane:'
read (*,*,err=7) vers1x, vers1y, vers1z
!_________________ if the plane is chosen by three points
else if(icase.eq.2) then
21 write(*,*)'Give the 1st (A) point:'
call givepoint(Ra(1),Ra(2),Ra(3),angstr)
write(*,*)'Give the 2nd (B) point:'
call givepoint(Rb(1),Rb(2),Rb(3),angstr)
Rbb(1)=Rb(1)-Ra(1)
Rbb(2)=Rb(2)-Ra(2)
Rbb(3)=Rb(3)-Ra(3)
write(*,*)'Give the 3rd (C) point:'
call givepoint(Rc(1),Rc(2),Rc(3),angstr)
Rcc(1)=Rc(1)-Ra(1)
Rcc(2)=Rc(2)-Ra(2)
Rcc(3)=Rc(3)-Ra(3)
vers1x=Rbb(2)*Rcc(3)-Rbb(3)*Rcc(2)
vers1y=Rbb(3)*Rcc(1)-Rbb(1)*Rcc(3)
vers1z=Rbb(1)*Rcc(2)-Rcc(1)*Rbb(2)
a=vers1x**2+vers1y**2+vers1z**2
if(a.lt.tiny) then
write(*,*)'Error! Your points lie on a line!'
write(*,*)'Return/Try again (r,R/other char)?'
read(*,'(a)') answer
if(answer.ne.'r' .and. answer.ne.'R') go to 21
return
end if
end if
!
!__________ normalize the normal and choose (rather arbitrarily)
! two perpendicular vectors in the plane
!
call normalize(vers1x,vers1y,vers1z)
if (vers1x.eq.0.and.vers1y.eq.0) then
vers2x=0.0
vers2y=-1.0
vers2z=0.0
else
vers2x=vers1y
vers2y=-vers1x
vers2z=0.0
end if
call normalize(vers2x,vers2y,vers2z)
if (vers1x.eq.0.and.vers1y.eq.0) then
vers3x=1.0
vers3y=0.0
vers3z=0.0
else
vers3x=-vers1x*vers1z
vers3y=-vers1y*vers1z
vers3z=vers1x**2+vers1y**2
end if
call normalize(vers3x,vers3y,vers3z)
end subroutine vector3
subroutine centralP(DIRC)
!.....................................................................
! Asks for the central point (acentx,acenty,acentz) in the plane
!.....................................................................
! icase = 1 - by a normal vector;
! 2 - by 3 points lying in the plane.
! angstr - gives a method of specifying coordinates
!.....................................................................
use menu
implicit none
real*8 DIRC(3,3),alpha
!
!............. give the central point of the plain: ................
!
!______ if icase=2 (the plane was given by 3 points) and icase1=2,
! the center of the triangle is chosen as the center of the plane
!
!______ in other cases the center is chosen in a general way and
! (if icase=2 and icase1=1) the projection of it is made
! to the plane with respect to the first point A (Ra):
!
write(*,*)'Give CENTRAL POINT on the plane:'
call givepoint(aCENTX,aCENTY,aCENTZ,angstr)
!
!________ removing a part normal to the plane (i.e. the projection)
if(icase.eq.2) then
write(*,*)'WARNING: The projection of this point on the'
write(*,*)' plane with respect to your A-point, '
write(*,'(a3,2(f10.5,a1),f10.5,a2)') &
' (',Ra(1),',',Ra(2),',',Ra(3),'),'
write(*,*)' will be taken as the central point.'
alpha=(acentx-Ra(1))*vers1x + (acenty-Ra(2))*vers1y + &
(acentz-Ra(3))*vers1z
acentx=acentx-alpha*vers1x
acenty=acenty-alpha*vers1y
acentz=acentz-alpha*vers1z
write(*,*)'Your central point (in Angstroms) is given by:'
write(*,'(a3,2(f10.5,a1),f10.5,a2)') &
' (',acentx,',',acenty,',',acentz,').'
end if
end subroutine centralP
subroutine givepoint2(aCENTX,aCENTY,angstr)
!...............................................................
! Choose a point for the plot on the plane
!...............................................................
use param
use atoms
implicit none
real*8 x(2),aCENTX,aCENTY
character angstr*12
integer number
!
!__________ give coordinates in Angstroms
if(angstr.eq.'<Angstroms> ') then
5 write(*,*)'Give Cartesian coordinates as X,Y in Angstroms:'
READ(*,*,err=5) aCENTX,aCENTY
!
!__________ give position via atomic number
else if(angstr.eq.'<AtomNumber>') then
30 write(*,*)'Give atomic number:'
read(*,*,err=30) number
if(number.lt.1.or.number.gt.NIONS) go to 31
aCENTX=TI(1,number)
aCENTY=TI(2,number)
return
31 write(*,*)'Error! Wrong atomic number! Try again!'
go to 30
!
!__________ give through fractional coordinates
else
6 write(*,*)'Give fractional coordinates along A1,A2:'
READ(*,*,err=11) x(1),x(2)
go to 20
11 write(*,*)'Error! Try again!'
go to 6
!_________ transform to Cartesian coordinates
20 aCENTX=x(1)*DIRC(1,1)+x(2)*DIRC(2,1)
aCENTY=x(1)*DIRC(1,2)+x(2)*DIRC(2,2)
return
end if
end subroutine givepoint2
subroutine givepoint(aCENTX,aCENTY,aCENTZ,angstr)
!...............................................................
! Choose a point for the plot
!...............................................................
use param
use atoms
implicit none
real*8 x(3),aCENTX,aCENTY,aCENTZ
character angstr*12
integer number
!
!__________ give coordinates in Angstroms
if(angstr.eq.'<Angstroms> ') then
5 write(*,*)'Give Cartesian coordinates as X,Y,Z in Angstroms:'
READ(*,*,err=5) aCENTX,aCENTY,aCENTZ
!
!__________ give position via atomic number
else if(angstr.eq.'<AtomNumber>') then
30 write(*,*)'Give atomic number:'
read(*,*,err=30) number
if(number.lt.1.or.number.gt.NIONS) go to 31
aCENTX=TI(1,number)
aCENTY=TI(2,number)
aCENTZ=TI(3,number)
return
31 write(*,*)'Error! Wrong atomic number! Try again!'
go to 30
!
!__________ give through fractional coordinates
else
6 write(*,*)'Give fractional coordinates along A1,A2,A3:'
READ(*,*,err=11) x(1),x(2),x(3)
go to 20
11 write(*,*)'Error! Try again!'
go to 6
!_________ transform to Cartesian coordinates
20 aCENTX=x(1)*DIRC(1,1)+x(2)*DIRC(2,1)+x(3)*DIRC(3,1)
aCENTY=x(1)*DIRC(1,2)+x(2)*DIRC(2,2)+x(3)*DIRC(3,2)
aCENTZ=x(1)*DIRC(1,3)+x(2)*DIRC(2,3)+x(3)*DIRC(3,3)
end if
end subroutine givepoint