← Back to branch summary

~nickpapior/siesta/trunk-kpoint-dos

« back to all changes in this revision

Viewing changes to Util/Denchar/Src/readpla.f

  • Committer: Alberto Garcia
  • Date: 2004-11-25 18:49:43 UTC
  • Revision ID: Arch-1:siesta@uam.es--2004%siesta-devel--reference--0.11--patch-1
Siesta 0.11 -- imported from CVS
Import from cvs using date instead of siesta-0-11-release tag, since
the Pseudo structure was not properly integrated at that time and
did not get the tag.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Util/Denchar/Src/readpla.f
 
1
 
 
2
      SUBROUTINE READPLA( IOPTION, XMIN, XMAX, YMIN, YMAX,
 
3
     .                    NPX, NPY, COORPO, NORMAL, DIRVER1, DIRVER2, 
 
4
     .                    ARMUNI, MAXA, XA, VOLUME, IUNITCD ) 
 
5
C **********************************************************************
 
6
C Read the data file to prepare the plane in which we are going to
 
7
C calculate the charge density
 
8
C Coded by J. Junquera, November 98
 
9
C **********************************************************************
 
10
 
 
11
      IMPLICIT NONE
 
12
 
 
13
      INCLUDE 'fdfdefs.h'
 
14
 
 
15
      INTEGER 
 
16
     .  IOPTION, NPX, NPY
 
17
     
 
18
      DOUBLE PRECISION
 
19
     .  XMIN, XMAX, YMIN, YMAX, 
 
20
     .  COORPO(3,3), NORMAL(3), DIRVER1(3), DIRVER2(3),
 
21
     .  ARMUNI
 
22
 
 
23
      INTEGER
 
24
     .  MAXA
 
25
      
 
26
      DOUBLE PRECISION
 
27
     .  XA(3,MAXA), VOLUME
 
28
 
 
29
C **********************************************************************
 
30
C INTEGER IOPTION        : Option to generate the plane
 
31
C                          1 = Normal vector
 
32
C                          2 = Two vectors belonging to the plane
 
33
C                          3 = Three points of the plane
 
34
C                          4 = Three atomic indices
 
35
C INTEGER NPX, NPY       : Number of points generated along x and y
 
36
C                          directions ina a system of reference in which
 
37
C                          the third component of the points of the plane
 
38
C                          is zero (Plane Reference Frame; PRF)
 
39
C REAL*8  XMIN, XMAX     : Limits of the plane in the PRF for x-direction
 
40
C REAL*8  YMIN, YMAX     : Limits of the plane in the PRF for y-direction
 
41
C REAL*8  NORMAL(3)      : Components of the normal vector used to define
 
42
C                          the plane 
 
43
C REAL*8  DIRVER(3)      : Components of the two vectors contained in the plane 
 
44
C REAL*8  COORPO(3,3)    : Coordinates of the three points used to define 
 
45
C                          the plane. COORPO(POINT,IX)
 
46
C INTEGER MAXA           : Maximum number of atoms
 
47
C REAL*8  XA(3,MAXA)     : Atomic coordinates
 
48
C REAL*8 VOLUME          : Volumen of unit cell (in bohr**3)
 
49
C **********************************************************************
 
50
 
 
51
      CHARACTER 
 
52
     .  OGP*22, OGP_DEFECT*22,
 
53
     .  CPF*22, CPF_DEFECT*22,
 
54
     .  UCD*22, UCD_DEFECT*22
 
55
 
 
56
      INTEGER
 
57
     .  IUNIT, IX, JX, NPX_DEFECT, NPY_DEFECT, IND1, IND2, IND3,
 
58
     .  ISCALE, IUNITCD
 
59
 
 
60
      DOUBLE PRECISION
 
61
     .  ORIGIN(3), XDIR(3)
 
62
 
 
63
      LOGICAL 
 
64
     .  LEQI, COLIN
 
65
 
 
66
      EXTERNAL 
 
67
     .  LEQI, COLINEAR
 
68
 
 
69
      DATA ORIGIN /0.D0,0.D0,0.D0/
 
70
      DATA XDIR   /1.D0,0.D0,0.D0/
 
71
      DATA IND1   /1/
 
72
      DATA IND2   /2/
 
73
      DATA IND3   /3/
 
74
      DATA COLIN  /.FALSE./
 
75
 
 
76
      CPF_DEFECT = 'Bohr'
 
77
      CPF = FDF_STRING('2D.CoorUnits',CPF_DEFECT)
 
78
      IF (LEQI(CPF,'bohr')) then
 
79
        ISCALE = 1
 
80
      ELSEIF (LEQI(CPF,'ang')) then
 
81
        ISCALE = 2
 
82
      ENDIF
 
83
 
 
84
      UCD_DEFECT = 'Ele/bohr**3'
 
85
      UCD = FDF_STRING('2D.DensityUnits',UCD_DEFECT)
 
86
      IF (LEQI(UCD,'ele/bohr**3')) then
 
87
        IUNITCD = 1
 
88
      ELSEIF (LEQI(UCD,'ele/ang**3')) then
 
89
        IUNITCD = 2
 
90
      ELSEIF (LEQI(UCD,'ele/unitcell')) then
 
91
        IUNITCD = 3
 
92
      ELSE
 
93
       WRITE(6,'(A)')' readpla:  Wrong Option in Units of             '
 
94
       WRITE(6,'(A)')' readpla:        Charge Density                 '
 
95
       WRITE(6,'(A)')' readpla:  You must choose one of the following:'       
 
96
       WRITE(6,'(A)')'                                                '
 
97
       WRITE(6,'(A)')' readpla:      - Ele/bohr**3                    '
 
98
       WRITE(6,'(A)')' readpla:      - Ele/ang**3                     '
 
99
       WRITE(6,'(A)')' readpla:      - Ele/unitcell                   '
 
100
       STOP
 
101
      ENDIF
 
102
 
 
103
      NPX_DEFECT = 50
 
104
      NPY_DEFECT = 50
 
105
      NPX = FDF_INTEGER('2D.NumberPointsX',NPX_DEFECT)
 
106
      NPY = FDF_INTEGER('2D.NumberPointsY',NPY_DEFECT)
 
107
 
 
108
      XMIN = FDF_PHYSICAL('2D.MinX',-3.D0,'Bohr')
 
109
      XMAX = FDF_PHYSICAL('2D.MaxX', 3.D0,'Bohr')
 
110
      YMIN = FDF_PHYSICAL('2D.MinY',-3.D0,'Bohr')
 
111
      YMAX = FDF_PHYSICAL('2D.MaxY', 3.D0,'Bohr')
 
112
 
 
113
      OGP_DEFECT = 'NormalVector'
 
114
      OGP = FDF_STRING('2D.PlaneGeneration',OGP_DEFECT)
 
115
      IF (LEQI(OGP,'normalvector')) then
 
116
        IOPTION = 1
 
117
      ELSEIF (LEQI(OGP,'twolines')) then
 
118
        IOPTION = 2
 
119
      ELSEIF (LEQI(OGP,'threepoints')) then
 
120
        IOPTION = 3
 
121
      ELSEIF (LEQI(OGP,'threeatomicindices')) then
 
122
        IOPTION = 4
 
123
      ELSE
 
124
       WRITE(6,'(A)')' readpla:  Wrong Option to Generate The Plane   '
 
125
       WRITE(6,'(A)')' readpla:  You must choose one of the following:'       
 
126
       WRITE(6,'(A)')'                                                '
 
127
       WRITE(6,'(A)')' readpla:      - NormalVector                   '
 
128
       WRITE(6,'(A)')' readpla:      - TwoLines                       '
 
129
       WRITE(6,'(A)')' readpla:      - ThreePoints                    '
 
130
       WRITE(6,'(A)')' readpla:      - ThreeAtomicIndices             '
 
131
       STOP
 
132
      ENDIF
 
133
 
 
134
      IF ( FDF_BLOCK('2D.CompNormalVector',IUNIT) ) THEN
 
135
        READ(IUNIT,*)(NORMAL(IX),IX=1,3)
 
136
      ENDIF
 
137
 
 
138
      IF ( FDF_BLOCK('2D.Comp2Vectors',IUNIT) ) THEN
 
139
        READ(IUNIT,*)(DIRVER1(IX),IX=1,3)
 
140
        READ(IUNIT,*)(DIRVER2(IX),IX=1,3)
 
141
      ENDIF
 
142
 
 
143
      IF ( FDF_BLOCK('2D.Coor3Points',IUNIT) ) THEN
 
144
        READ(IUNIT,*)(COORPO(1,IX),IX=1,3)
 
145
        READ(IUNIT,*)(COORPO(2,IX),IX=1,3)
 
146
        READ(IUNIT,*)(COORPO(3,IX),IX=1,3)
 
147
      ENDIF
 
148
 
 
149
      IF ( FDF_BLOCK('2D.Indices3Atoms',IUNIT) ) THEN
 
150
        READ(IUNIT,*)IND1, IND2, IND3
 
151
      ENDIF
 
152
 
 
153
      IF ( IOPTION .EQ. 4 ) THEN
 
154
        DO IX = 1,3
 
155
          COORPO(1,IX) = XA(IX,IND1)
 
156
        ENDDO
 
157
        DO IX = 1,3
 
158
          COORPO(2,IX) = XA(IX,IND2)
 
159
        ENDDO
 
160
        DO IX = 1,3
 
161
          COORPO(3,IX) = XA(IX,IND3)
 
162
        ENDDO
 
163
      ENDIF
 
164
 
 
165
C Check if the three points are colinear -------------------------------
 
166
      IF ((IOPTION .EQ. 3) .OR. (IOPTION .EQ. 4))THEN
 
167
         CALL COLINEAR( COORPO, COLIN )
 
168
         IF(COLIN) THEN
 
169
           WRITE(6,*)'The coordinates of the three points are colinear'
 
170
           WRITE(6,*)'and do not define a plane' 
 
171
           WRITE(6,*)'Please, check these coordinates in the input file'      
 
172
           STOP
 
173
         ENDIF
 
174
      ENDIF
 
175
 
 
176
      IF ( FDF_BLOCK('2D.PlaneOrigin',IUNIT) ) THEN
 
177
        READ(IUNIT,*)(ORIGIN(IX),IX=1,3)
 
178
      ENDIF
 
179
 
 
180
      IF ( FDF_BLOCK('2D.X_Axis',IUNIT) ) THEN
 
181
        READ(IUNIT,*)(XDIR(IX),IX=1,3)
 
182
      ENDIF
 
183
 
 
184
      IF (IOPTION .LT. 3) THEN
 
185
        DO IX = 1,3      
 
186
          COORPO(1,IX) = ORIGIN(IX)
 
187
        ENDDO
 
188
        IF(IOPTION .EQ. 1) THEN
 
189
          DO IX = 1,3
 
190
            COORPO(2,IX) = XDIR(IX)
 
191
          ENDDO
 
192
        ENDIF
 
193
      ENDIF
 
194
 
 
195
C Scale points coordinates
 
196
C   Iscale = 1 => Do nothing
 
197
C   Iscale = 2 => Multiply by 1./0.529177 (Ang --> Bohr)
 
198
 
 
199
      IF( (ISCALE .EQ. 2) .AND. (IOPTION .NE. 4) ) THEN
 
200
        DO IX = 1,3
 
201
          DO JX = 1,3
 
202
            COORPO(JX,IX) = 1.D0 / 0.529177D0 * COORPO(JX,IX)
 
203
          ENDDO
 
204
          ORIGIN(IX)  = 1.D0 / 0.529177D0 * ORIGIN(IX)
 
205
          DIRVER1(IX) = 1.D0 / 0.529177D0 * DIRVER1(IX)
 
206
          DIRVER2(IX) = 1.D0 / 0.529177D0 * DIRVER2(IX)
 
207
        ENDDO
 
208
      ENDIF 
 
209
 
 
210
C Units of Charge Density
 
211
C   Iunitcd = 1 => Do nothing
 
212
C   Iunitcd = 2 => Multiply by (1.d0 / 0.529177d0) **3 (bohr**3 --> Ang**3)
 
213
C   Iunitcd = 3 => Multiply by volume unit cell (in bohrs**3) 
 
214
 
 
215
      IF (IUNITCD .EQ. 1) THEN
 
216
        ARMUNI = 1.D0
 
217
      ELSEIF( IUNITCD .EQ. 2 ) THEN
 
218
        ARMUNI = (1.D0 / 0.529177D0)**3 
 
219
      ELSEIF( IUNITCD .EQ. 3 ) THEN
 
220
        ARMUNI = VOLUME
 
221
      ENDIF
 
222
 
 
223
      END
 
224