1 PROGRAM GENGEOM 2C ******************************************* 3C *** GENERAL GEOMETRIC PURPOSE PROGRAM *** 4C ******************************************* 5 6c ======================================================================= 7c USAGE: 8c gengeom MODE1 MODE2 MODE3 IGRP NXDIR NYDIR NZDIR OUTPUT INPUT 9c 0 1 2 3 4 5 6 7 8 9 10c ======================================================================= 11 12c program read data from ftn34, but if INPUT ARGUMENT is present it 13c rather read data from gengeom_formated INPUT file instead 14 15C NXDIR,NYDIR,NZDIR --> n. of repeatitions in each directions 16 17C ***** 18C IGR - SEQUENTIAL GROUP NUMBER (USED ONLY FOR CRYSTALS) 19C posebej je treba paziti na R/H !!!!! 20C ***** 21 22c how the vectors are packed to matrix:: 23c -------------------------------------- 24c DIRECT VECTORS:: 25c ---------------- 26c VECTOR1 GOES TO FIRST COLOUMN, VETCOR2 TO SECOND & VECTOR3 TO 27c THIRD COLOUMN;;; DV(COL,ROW); vec(# ,x/y/z) 28c ============= 29c |11 21 31| |x1 x2 x3| 30c VEC = |12 22 32| => DV = |y1 y2 y3| 31c |13 23 33| |z1 z2 z3| 32c 33c RECIPROCAL VECTORS:: (in transpose manner) 34c -------------------- 35c VECTOR1 GOES TO FIRST ROW, VETCOR2 TO SECOND & VECTOR3 TO 36c THIRD ROW;;; vec(x/y/z, #) 37c ============= 38c |11 21 31| |x1* y1* z1*| 39c IVEC = |12 22 32| =>IDV = |x2* y2* z2*| 40c |13 23 33| |x3* y3* z3*| 41c 42c so: |1 0 0| 43c VEC x IVEC == |0 1 0| 44c |0 0 1| 45 46 IMPLICIT REAL*8 (A-H,O-Z) 47 include 'param.inc' 48 CHARACTER*256 MODE1, MODE2, MODE3, 49 * ANX, ANY, ANZ, FILE2, INPUT 50 REAL*8 PC(3,1),AC(3,2),BC(3,2),CC(3,2),FC(3,4),IC(3,2),RC(3,3), 51 * HC(3,3),CSTMC(3,4),P2A(3,3),P2B(3,3),P2C(3,3),P2F(3,3), 52 * P2I(3,3),R2H(3,3),P2AI(3,3),P2BI(3,3),P2CI(3,3),P2FI(3,3), 53 * P2II(3,3),R2HI(3,3),P2H(3,3),P2HI(3,3), 54 * DV(3,3),SOP(3,3,48),TRX(48),TRY(48),TRZ(48), 55 * DVC(3,3),CON13,CON23,HTTC(3,13),RTTC(3,13),MINTOL, 56 * DVEC(3,3),TMPVEC(3,3),X(NAC),Y(NAC),Z(NAC), 57 * IDVEC(3,3),IDV(3,3),IDVC(3,3), 58 * RCR(3,3),RTTCR(3,13),S_AC(3,2),S_BC(3,2),S_CC(3,2), 59 * S_FC(3,4),S_IC(3,2),S_RC(3,3),S_HC(3,3) 60 INTEGER C2I, NAT(NAC) 61 LOGICAL L, lreverse, IsReverse, notfound, Lprimvec 62 63 include 'mode.inc' 64 65 COMMON/FTN34/ SOP,TRX,TRY,TRZ,DV,IDIM,IGROUP,NSYMOP 66 COMMON/APOS/ AC,BC,CC,FC,IC,HC,RC,DVC,CSTMC 67 COMMON/DIR/ NXDIR,NYDIR,NZDIR 68 COMMON/IVEC/ IDVEC 69 COMMON/MULTAT/ X,Y,Z,NAT,DVEC,NATR 70 COMMON/BOHRR/ BOHR,IMODE3 71 COMMON/KEYWORD/ notfound, Lprimvec 72 73 PARAMETER (CON13=1.0d0/3.0d0, CON23=2.0d0/3.0d0, MINTOL=1.d-6, 74 $ CON13N=-1.0d0/3.0d0, CON23N=-2.0d0/3.0d0, 75 $ TOLMIN=1.d-6) 76 77C VARIABLE NAMES: 78C DV{1|2|3}.............3 direct vectors in PRIMITIV cell 79C DVC{1|2|3}............3 direct vectors in CONVENTIONAL cell 80C TR{X|Y|Z}.............X,Y,Z COMPONENT OF TRANSL. DUE TO POINT SYM. OPER. 81C IDIM..................DIMENSION OF THE SYSTEM 82C IGROUP................GROUP NUMBER (LOOK BELOW) 83C SOP...................SYM. OPERATORS (MATRICES) 84 85C *********************************************** 86C TRANSFORMATION of PRIMITIV CELL to CONVENTIONAL 87C *********************************************** 88C =================== 89C IGROUP DEFINITIONS: 90C =================== 91C GROUP NUMBER USED IN ftn34 92C N. symbol N.of atoms Variable 93C 1.....P --> 1 PC 94C 2.....A --> 2 AC 95C 3.....B --> 2 BC 96C 4.....C --> 2 CC 97C 5.....F --> 4 FC 98C 6.....I --> 2 IC 99C 7.....R --> 3 RC --> RTTC 100C REMARK: ALL HEXAGONAL CELLS ARE OF TYPE 1 IN FTN34 !!!!! 101C SO I'LL USE 8 FOR HEXAGONAL CELL 102C 8.....H --> 3 HC --> HTTC 103C for TRIGONAL CELL that are not RHOMOHEDRAL, I will use additional flag 104C 9.....TRIGONAL_NOT_RHOMOHEDRAL 105C 10....custom; used just when calling CONVATOM to tell to read CSTMC 106C ---------------------------------------- 107C ---------------------------------------- 108C TRANSFORMATION MATRICES: 109C ------------------------ 110C EXAMPLE: 111C P2C means: PRIMITIV to C-centered 112C P2C and the rest refer to coor. of point, so 113C P2CI and the rest refer to transf. of vectors 114C P2CI means "inverse of P2C" 115 116 DATA PC/ 117 1 0.0d0, 0.0d0, 0.0d0/ 118C CONVENTIONAL CELL; COORD. OF LATTICE POINTS WITHIN CELL 119 DATA S_AC/ 120 1 0.0d0, 0.0d0, 0.0d0, 121 2 0.0d0, 0.5d0, 0.5d0/ 122 DATA S_BC/ 123 1 0.0d0, 0.0d0, 0.0d0, 124 2 0.5d0, 0.0d0, 0.5d0/ 125 DATA S_CC/ 126 1 0.0d0, 0.0d0, 0.0d0, 127 2 0.5d0, 0.5d0, 0.0d0/ 128 DATA S_FC/ 129 1 0.0d0, 0.0d0, 0.0d0, 130 2 0.5d0, 0.5d0, 0.0d0, 131 3 0.0d0, 0.5d0, 0.5d0, 132 4 0.5d0, 0.0d0, 0.5d0/ 133 DATA S_IC/ 134 1 0.0d0, 0.0d0, 0.0d0, 135 2 0.5d0, 0.5d0, 0.5d0/ 136 DATA S_HC/ 137 1 0.0d0, 0.0d0, 0.0d0, 138 2 CON23, CON13, 0.0d0, 139 3 CON13, CON23, 0.0d0/ 140C RHOMBOHEDRALY CENTERED (OBVERSE SETTINGS) - hexagonal axes 141C CRYSTALXX uses this ^^^^^^^ 142 DATA S_RC/ 143 1 0.0d0, 0.0d0, 0.0d0, 144 2 CON23, CON13, CON13, 145 3 CON13, CON23, CON23/ 146 147C RHOMBOHEDRALY CENTERED (REVERSE SETTINGS) - hexagonal axes 148C WIEN uses this ^^^^^^^ 149 DATA RCR/ 150 1 0.0d0, 0.0d0, 0.0d0, 151 2 CON13, CON23, CON13, 152 3 CON23, CON13, CON23/ 153 154C HEXAGONALY TRIPLE CENTERED CELL - HEXAGONAL SHAPE 155 DATA HTTC/ 156 1 1.0d0, 0.0d0, 0.0d0, 1.0d0, 1.0d0, 0.0d0, 157 2 0.0d0, 1.0d0, 0.0d0, -1.0d0, 0.0d0, 0.0d0, 158 3 -1.0d0,-1.0d0, 0.0d0, 0.0d0,-1.0d0, 0.0d0, 159 4 0.0d0, 0.0d0, 0.0d0, 160 5 CON23, CON13, 0.0d0, CON13, CON23, 0.0d0, 161 6 CON13N, CON13, 0.0d0, CON23N,CON13N, 0.0d0, 162 7 CON13N,CON23N, 0.0d0, CON13,CON13N, 0.0d0/ 163C RHOMBOHEDRAL - described as TRIPLE HEXAGONAL CELL - HEXAGONAL SHAPE 164c OBVERSE SETTINGS 165c ^^^^^^^ 166 DATA RTTC/ 167 1 1.0d0, 0.0d0, 0.0d0, 1.0d0, 1.0d0, 0.0d0, 168 2 0.0d0, 1.0d0, 0.0d0, -1.0d0, 0.0d0, 0.0d0, 169 3 -1.0d0,-1.0d0, 0.0d0, 0.0d0,-1.0d0, 0.0d0, 170 4 0.0d0, 0.0d0, 0.0d0, 171 5 CON23, CON13, CON13, CON13, CON23, CON23, 172 6 CON13N, CON13, CON13, CON23N,CON13N, CON23, 173 7 CON13N,CON23N, CON13, CON13,CON13N, CON23/ 174c REVERSE SETTINGS 175c ^^^^^^^ 176 DATA RTTCR/ 177 1 1.0d0, 0.0d0, 0.0d0, 1.0d0, 1.0d0, 0.0d0, 178 2 0.0d0, 1.0d0, 0.0d0, -1.0d0, 0.0d0, 0.0d0, 179 3 -1.0d0,-1.0d0, 0.0d0, 0.0d0,-1.0d0, 0.0d0, 180 4 0.0d0, 0.0d0, 0.0d0, 181 5 CON13, CON23, CON13, CON23, CON13, CON23, 182 6 CON23N, CON23, CON13, CON13N,CON23N, CON23, 183 7 CON23N,CON13N, CON13, CON23,CON23N, CON23/ 184 185C TRANSFORMATION MATRICES IN CRYSTAL95 186 DATA P2A/ 187 1 1.0d0, 0.0d0, 0.0d0, 188 2 0.0d0, 0.5d0,-0.5d0, 189 3 0.0d0, 0.5d0, 0.5d0/ 190 DATA P2AI/ 191 1 1.0d0, 0.0d0, 0.0d0, 192 2 0.0d0, 1.0d0, 1.0d0, 193 3 0.0d0,-1.0d0, 1.0d0/ 194 195 DATA P2B/ 196 1 0.5d0, 0.0d0, 0.5d0, 197 2 0.0d0, 1.0d0, 0.0d0, 198 3 -0.5d0, 0.0d0, 0.5d0/ 199 DATA P2BI/ 200 1 1.0d0, 0.0d0,-1.0d0, 201 2 0.0d0, 1.0d0, 0.0d0, 202 3 1.0d0, 0.0d0, 1.0d0/ 203 204 DATA P2C/ 205 1 0.5d0,-0.5d0, 0.0d0, 206 2 0.5d0, 0.5d0, 0.0d0, 207 3 0.0d0, 0.0d0, 1.0d0/ 208 DATA P2CI/ 209 1 1.0d0, 1.0d0, 0.0d0, 210 2 -1.0d0, 1.0d0, 0.0d0, 211 3 0.0d0, 0.0d0, 1.0d0/ 212 213 DATA P2F/ 214 1 0.0d0, 0.5d0, 0.5d0, 215 2 0.5d0, 0.0d0, 0.5d0, 216 3 0.5d0, 0.5d0, 0.0d0/ 217 DATA P2FI/ 218 1 -1.0d0, 1.0d0, 1.0d0, 219 2 1.0d0,-1.0d0, 1.0d0, 220 3 1.0d0, 1.0d0,-1.0d0/ 221 222 DATA P2I/ 223 1 -0.5d0, 0.5d0, 0.5d0, 224 2 0.5d0,-0.5d0, 0.5d0, 225 3 0.5d0, 0.5d0,-0.5d0/ 226 DATA P2II/ 227 1 0.0d0, 1.0d0, 1.0d0, 228 2 1.0d0, 0.0d0, 1.0d0, 229 3 1.0d0, 1.0d0, 0.0d0/ 230 231 DATA R2H/ 232 1 CON23,CON13N,CON13N, 233 2 CON13, CON13,CON23N, 234 3 CON13, CON13, CON13/ 235 DATA R2HI/ 236 1 1.0d0, 0.0d0, 1.0d0, 237 2 -1.0d0, 1.0d0, 1.0d0, 238 3 0.0d0,-1.0d0, 1.0d0/ 239C HEXAGONAL CELL P --> TRIPLE HEXAGONAL CELL H1; 240C INTERNATIONAL TABLES OF CRYSTALOGRAPHY, 1993; TABLE 5-1 241 DATA P2H/ 242 1 CON23,CON13N, 0.0d0, 243 2 CON13, CON13, 0.0d0, 244 3 0.0d0, 0.0d0, 1.0d0/ 245 DATA P2HI/ 246 1 1.0d0, 1.0d0, 0.0d0, 247 2 -1.0d0, 2.0d0, 0.0d0, 248 3 0.0d0, 0.0d0, 1.0d0/ 249C === END_OF_DATA_BLOCK === 250 251 call MATCOPY(S_AC,AC,3,2) 252 call MATCOPY(S_BC,BC,3,2) 253 call MATCOPY(S_CC,CC,3,2) 254 call MATCOPY(S_FC,FC,3,4) 255 call MATCOPY(S_IC,IC,3,2) 256 call MATCOPY(S_RC,RC,3,3) 257 call MATCOPY(S_HC,HC,3,3) 258 259C CONVERSION FROM BOHRS TO ANGS. 260 BOHR=0.529177d0 261 262 NCSTMC=0 263 264c number of command line arguments must not be lower than 8 265 NARG=IARGC() 266 IF(NARG.LT.8)THEN 267 WRITE(*,*) 268 1 'usage: gengeom <mode1> <mode2> <mode3> <igrp> <nxdir> <n 269 2ydir> <nzdir> <output> [input]' 270 STOP 271 ENDIF 272 273C *********************************** 274C *** READ MODE & N?DIR ARGUMENTS *** 275C *********************************** 276C look in 'mode.inc' for definitions !!!! 277 CALL GETARG(IARG_MODE1,MODE1) 278 CALL GETARG(IARG_MODE2,MODE2) 279 CALL GETARG(IARG_MODE3,MODE3) 280 CALL GETARG(IARG_NXDIR,ANX) 281 CALL GETARG(IARG_NYDIR,ANY) 282 CALL GETARG(IARG_NZDIR,ANZ) 283 IMODE1=C2I(MODE1) 284 IMODE2=C2I(MODE2) 285 IMODE3=C2I(MODE3) 286 NXDIR =C2I(ANX) 287 NYDIR =C2I(ANY) 288 NZDIR =C2I(ANZ) 289c determine aproximate number of atoms (COUNT JUST ONE ATOM PER CELL) 290 NATOMS=(NXDIR+2) * (NYDIR+2) *(NZDIR+1) !overestimated, thats good 291 292 CALL GETARG(IARG_OUTPUT,FILE2) 293 IFILE2=INDEX(FILE2,' ') 294 OPEN(UNIT=12,FILE=FILE2(1:IFILE2-1),STATUS='UNKNOWN') 295C file2='tmp.tmp' 296C OPEN(UNIT=12,FILE='tmp.tmp',STATUS='UNKNOWN') 297 298c --- DEBUG_BEGIN --- 299c print *,'IMODE1=',IMODE1,'; IMODE2=',IMODE2,'; IMODE3=',IMODE3 300c print *,'NXDIR=',NXDIR,'; NYDIR=',NYDIR,'; NZDIR=',NZDIR 301c print *,'IFILE2=',IFILE2,'; OUTPUT::',FILE2(1:IFILE2-1) 302c --- DEBUG_END --- 303 304c ********************************* 305c take care of imode3; imode3 comprise several modes, that are based on 306c following criteria: 1258, 1 is mode3, 2 is mode4, 5 is mode5, 8 is mode6 307c so far only MODE3 & MODE4 are used!!!! 308 IMODE4=MOD(IMODE3,10) 309 IMODE3=(IMODE3-IMODE4)/10 310 311C ********************************* 312C *** write INFO record *** 313 WRITE(12,'(1x,a)') 'INFO' 314 WRITE(12,'(a,3x,i4,1x,i4,1x,i4)') 'nunit',NXDIR,NYDIR,NZDIR 315 IF(IMODE2.EQ.M2_CELL) WRITE(12,'(a,3x,a)') 'unit','cell' 316 IF(IMODE2.EQ.M2_TR_ASYM_UNIT) 317 1 WRITE(12,'(a,3x,a)') 'unit','tr_asym' 318 IF(IMODE1.EQ.M1_PRIM .OR. IMODE1.EQ.M1_PRIM3) THEN 319 WRITE(12,'(a,3x,a)') 'celltype', 'primcell' 320 ELSE 321 WRITE(12,'(a,3x,a)') 'celltype', 'convcell' 322 ENDIF 323 IF(IMODE1.EQ.M1_HEXA_SHAPE .OR. IMODE1.EQ.M1_PRIM3) THEN 324 WRITE(12,'(a,3x,a)') 'shape','hexagonal' 325 ELSE 326 WRITE(12,'(a,3x,a)') 'shape','parapipedal' 327 ENDIF 328 WRITE(12,'(1x,a)') 'END_INFO' 329C ********************************* 330 331 IGROUP_input=0 332 IF(NARG.EQ.9)THEN 333C ********************************* 334C *** is INPUT ARGUMENT PRESENT *** 335C *** read from INPUT *** 336C ********************************* 337 IDUM1=0 338 IDUM2=0 339 CALL GETARG(IARG_INPUT,INPUT) 340 I_INPUT=INDEX(INPUT,' ') 341 OPEN(UNIT=11,FILE=INPUT(1:I_INPUT),STATUS='OLD') 342c READF1 has two common's that are: 343c /MULTAT/ X,Y,Z,NAT,DVEC,NATR 344c /IVEC/ IDVEC 345 346c backward compatibility 347 CALL READF1('DIM-GROUP',IGROUP,IDIM) 348c -- 349 CALL READF1('MOLECULE',IGROUP,IDIM) 350 CALL READF1('POLYMER',IGROUP,IDIM) 351 CALL READF1('SLAB',IGROUP,IDIM) 352 CALL READF1('CRYSTAL',IGROUP,IDIM) 353 354 CALL READF1('PRIMVEC',IGROUP,IDIM) 355 if (notfound) then 356 Lprimvec=.false. 357 else 358 Lprimvec=.true. 359 CALL MATCOPY(DVEC,DV,3,3) 360 endif 361 362 CALL READF1('CONVVEC',IGROUP,IDIM) 363 if ((.not.notfound) .and. (.not.Lprimvec) ) then 364 call Matcopy(dvec,dv,3,3) 365 write(12,'(1x,a)') 'PRIMVEC' 366 write(12,'(3(1x,f15.10)/3(1x,f15.10)/3(1x,f15.10))') 367 1 ((dvec(j,i),i=1,3),j=1,3) 368 else if (notfound .and. (.not.Lprimvec)) then 369 WRITE(6,*) 'PRIMVEC and CONVEC keywords not found !!!' 370 STOP 371 endif 372 CALL MATCOPY(DVEC,DVC,3,3) 373 374 if (Lprimvec) then 375 CALL READF1('PRIMCOORD',IDUM1,IDUM2) !here NAT is read 376 else 377c assuming PRIMCOORD == CONVCOORD 378 CALL READF1('CONVCOORD',IDUM1, IDUM2) 379 endif 380 381 CALL GetCenteredPoints(DV, DVC, IDIM, CSTMC, NCSTMC) 382 IGROUP_input=IGROUP !to save original igroup read from input 383 IGROUP=10 !10 means custom (input file was read) 384 385c DEBUG_BEGIN 386c if (IGROUP_input.EQ.7) 387c $ print *,'DEBUG: IsReverse = ', IsReverse(CSTMC,RCR) 388c print *,'DEBUG: NCSTMC::', NCSTMC 389c do i=1,4 390c print *,'DEBUG: CSTMC::',(CSTMC(j,i),j=1,3) 391c enddo 392c DEBUG_END 393 394 if (ncstmc.gt.4) then 395c Largest number of special centered points has FCC == 4; if NCSTMC exceed 396c 4, something is wrong 397 STOP 'more then four centered points; conventional vectors 398 $ are probably badly chosen' 399 ENDIF 400 CLOSE(11) 401 ELSE 402C ************************************************* 403C ***** READ FTN34; all modes need this ******* 404C ************************************************* 405 NATR=0 406c READFTN34 has one common: 407C /FTN34/ IDIM,IGROUP,DV,NSYMOP,SOP,TRX,TRY,TRZ 408 CALL READFTN34(NATR, IMODE4) !NATR is read 409 410C **************************** 411C get the CONVENTIONAL VECTORS 412C **************************** 413 call ZEROMAT(DVC, 3, 3) 414C P->P 415 IF(IGROUP.EQ.1) THEN 416 DO J=1,3 417 DO I=1,3 418 DVC(I,J)=DV(I,J) 419 ENDDO 420 ENDDO 421 ENDIF 422C P->A 423 IF(IGROUP.EQ.2) CALL MATMULT(DV,3,3,P2AI,3,DVC) 424C P->B 425 IF(IGROUP.EQ.3) CALL MATMULT(DV,3,3,P2BI,3,DVC) 426C P->C 427 IF(IGROUP.EQ.4) CALL MATMULT(DV,3,3,P2CI,3,DVC) 428C P->F 429 IF(IGROUP.EQ.5) CALL MATMULT(DV,3,3,P2FI,3,DVC) 430C P->I 431 IF(IGROUP.EQ.6) CALL MATMULT(DV,3,3,P2II,3,DVC) 432C PR->H; primitiv-rhombohedrall to triple-hexagonal R cell 433 IF(IGROUP.EQ.7) CALL MATMULT(DV,3,3,R2HI,3,DVC) 434C PH->H; primitiv hexagonal --> triple hexagonal H cell 435 IF(IGROUP.EQ.8 .OR. IGROUP.EQ.9) 436 1 CALL MATMULT(DV,3,3,P2HI,3,DVC) 437 438C WRITE CONVENTIONAL DIRECT VECTORS 439 CALL WVEC('CONVVEC',DVC) 440 ENDIF 441 442c /////////////////////////////////////////////////////// 443c // below is common to ftn34 and XCrySDen format file // 444c /////////////////////////////////////////////////////// 445 446C NOW CALCULATE POSITIONS OF ATOMS NEEDED WITHIN CONV.CELL 447 CALL CONVATOM(IGROUP,NARG,NCSTMC,IMODE4) !primcoord & convcoord are writen to UNIT12 448 449C make reciprocal vectors (TRANSPOSE form) 450 CALL ZeroMat(IDV, 3, 3) 451 CALL ZeroMat(IDVC, 3, 3) 452c initialise tmpvec matrix; tmpvec is dummy anyway 453 call InvertSqMat123(DV, IDV, IDIM) 454 call InvertSqMat123(DVC, IDVC, IDIM) 455 456c write PRIMITIV & CONVENTIONAL RECIPROCAL VECTORS in normal not 457c transpose form 458 CALL WIVEC('RECIP-PRIMVEC',IDV) 459 CALL WIVEC('RECIP-CONVVEC',IDVC) 460 461c write WIGNER-SEITZ CELL 462 IF(IDIM.EQ.3 .AND. IMODE1.NE.M1_INFO) THEN 463 CALL WignerSeitz(DV, 'WIGNER-SEITZ-PRIMCELL') 464 CALL WignerSeitz(DVC, 'WIGNER-SEITZ-CONVCELL') 465 Call MatTranspose(IDV, TMPVEC, 3, 3) 466 CALL WignerSeitz(TMPVEC, 'BRILLOUIN-ZONE-PRIMCELL') 467 Call MatTranspose(IDVC, TMPVEC, 3, 3) 468 CALL WignerSeitz(TMPVEC, 'BRILLOUIN-ZONE-CONVCELL') 469 ELSEIF(IDIM.EQ.2 .AND. IMODE1.NE.M1_INFO) THEN 470 CALL WignerSeitz2D(DV, 'WIGNER-SEITZ-PRIMCELL') 471 Call MatTranspose(IDV, TMPVEC, 3, 3) 472 CALL WignerSeitz2D(TMPVEC, 'BRILLOUIN-ZONE-PRIMCELL') 473 ELSEIF(IDIM.EQ.1 .AND. IMODE1.NE.M1_INFO) THEN 474 CALL WignerSeitz1D(DV, 'WIGNER-SEITZ-PRIMCELL') 475 CALL WignerSeitz1D(IDV, 'BRILLOUIN-ZONE-PRIMCELL') 476 ENDIF 477 478C **************************** 479C deside upon IMODE1 480 IF(IMODE1.EQ.M1_INFO) GOTO 1111 !info mode 481 GOTO(150,250,250,250) IMODE1 482C ============================ 483 150 CONTINUE 484C ******************************************************* 485C * CALCULATE COORDINATES FOR PRIMITIV CELL(S); * 486C * ^^^^^^^^ * 487C * IMODE1 = M1_PRIM * 488C ******************************************************* 489 490c ASSING primitiv-vectors to DVEC & IDVEC 491 CALL MATCOPY(DV,DVEC,3,3) !dvec is used in MULTATOM; common/MULTAT/ 492 CALL MATCOPY(IDV,IDVEC,3,3) 493 494C IF IDIM=0 JUST PRINT THE COORDINATES 495 IF(IDIM.EQ.0)THEN 496 CALL MOLECULE 497 ELSE 498c print *,'NATOMS',NATOMS 499 CALL MULTATOM(PC,1,1,NATOMS,.false.,IDIM,IMODE2) 500 ENDIF 501 502 GOTO 1111 503 504 250 CONTINUE 505C ********************************************************* 506C ** calculate coordinates for CONVENTIONAL cell(s) && ** 507C ** also M1_PRIM3; look below ^^^^^^^^^^^^ ** 508C ********************************************************* 509c 510C ********************************************************** 511C for TRIGONAL/HEXAGONAL/RHOMBOHEDRAL SYSTEMS: 512C ** IMODE1=M1_CONV ........ PARAPIPEDAL SHAPED TRIPLE hexagonal CELL 513C ** IMODE1=M1_HEXA_SHAPE .. HEXAGONAL SHAPED TRIPLE hexagonal CELL 514c 515C for TRIGONAL_NOT_RHOMBO/HEXAGONAL systems 516C ** IMODE1=M1_PRIM3 --> HEXAGONAL SHAPED PRIMITIV CELL 517c ^^^^^^^^ 518C ********************************************************** 519 520C *** so far, only crystals can be converted from prim. to conv. cell *** 521 IF(IDIM.NE.3) 522 1 STOP 'conversion to conventional cell for non-crystal' 523 524c ASSING conventional-direct-vectors to DVEC 525 CALL MATCOPY(DVC,DVEC,3,3) !dvec is used in MULTATOM & MULTHEXA 526 CALL MATCOPY(IDVC,IDVEC,3,3) 527 528C HEXAGONAL SHAPED PRIMITIV RHOMOHEDRAL CELL 529c IF(IMODE1.EQ.M1_PRIM3 .AND. IGROUP.EQ.7) !maybe turn this off 530c 1 CALL MULTATOM(RC,3,8,NATOMS*3,.true.,3,IMODE2) 531 532C HEXAGONAL SHAPED PRIMITIV TRIGONAL_NOT_RHOMBO/HEXAGONAL CELL 533C t.k.: maybe TRIGONAL is not allowed for this 534 IF(IMODE1.EQ.M1_PRIM3 .AND. 535 $ (IGROUP.EQ.8 .OR. IGROUP.EQ.9 .OR. 536 $ IGROUP_input.EQ.8 .OR. IGROUP_input.EQ.9)) then 537 CALL MULTATOM(HC,3,8,NATOMS*3,.true.,3,IMODE2) 538 goto 1111 539 ENDIF 540 541C WHEATER (NOT HEXA/RHOMBO/TRIGONAL_NOT_RHOMBO) 542 L=.FALSE. 543 IF(IGROUP.EQ.1) CALL MULTATOM(PC,1,IGROUP,NATOMS,L,3,IMODE2) 544 IF(IGROUP.EQ.2) CALL MULTATOM(AC,2,IGROUP,NATOMS*2,L,3,IMODE2) 545 IF(IGROUP.EQ.3) CALL MULTATOM(BC,2,IGROUP,NATOMS*2,L,3,IMODE2) 546 IF(IGROUP.EQ.4) CALL MULTATOM(CC,2,IGROUP,NATOMS*2,L,3,IMODE2) 547 IF(IGROUP.EQ.5) CALL MULTATOM(FC,4,IGROUP,NATOMS*4,L,3,IMODE2) 548 IF(IGROUP.EQ.6) CALL MULTATOM(IC,2,IGROUP,NATOMS*2,L,3,IMODE2) 549 IF(IGROUP.EQ.10) then 550 if(IGROUP_input.EQ.7 .and. IMODE1.EQ.M1_HEXA_SHAPE) then 551c rhombohedral 552 lreverse=IsReverse(CSTMC, RCR) 553 if(lreverse) then 554c reverse settings 555 CALL MULTHEXA(RTTCR,IGROUP,NATOMS*13) 556 goto 1111 557 endif 558 elseif(IGROUP_input.GE.8 .and. IMODE1.EQ.M1_HEXA_SHAPE) then 559c hexagonal 560 IGROUP=IGROUP_input 561 else 562 CALL MULTATOM(CSTMC,NCSTMC,IGROUP,NATOMS*NCSTMC,L,3,IMODE2) 563 goto 1111 564 endif 565 ENDIF 566 567C TRIPLE HEXA/RHOMBO 568C M1_HEXA_SHAPE & M2_TR_ASYM_UNIT together are not allowed 569 IF(IMODE1.EQ.M1_HEXA_SHAPE .AND. IMODE2.EQ.M2_TR_ASYM_UNIT) 570 * STOP 'M1_HEXA_SHAPE & M2_TR_ASYM_UNIT can not be used together' 571 IF(IMODE1.EQ.M1_CONV .AND. IGROUP.EQ.7) 572 * CALL MULTATOM(RC,3,7,NATOMS*3,L,3,IMODE2) 573 IF(IMODE1.EQ.M1_CONV .AND. (IGROUP.EQ.8 .OR. IGROUP.EQ.9)) 574 * CALL MULTATOM(HC,3,8,NATOMS*3,L,3,IMODE2) 575 IF(IMODE1.EQ.M1_HEXA_SHAPE .AND. IGROUP.EQ.7) 576 * CALL MULTHEXA(RTTC,IGROUP,NATOMS*13) 577 IF(IMODE1.EQ.M1_HEXA_SHAPE .AND. (IGROUP.EQ.8 .OR. IGROUP.EQ.9)) 578 * CALL MULTHEXA(HTTC,IGROUP,NATOMS*13) 579 580 1111 CONTINUE 581 CLOSE(12) 582 END 583 584 585c INTEGER FUNCTION DETGROUP(IGR) 586C FUNCTION DETERMINES IF GROUP IS OF HEXAGONAL TYPE 587C 8 FOR HEXAGONAL TYPE 588c DETGROUP=1 589c IF(IGR.GE.143.AND.IGR.LE.194)DETGROUP=8 590c RETURN 591c END 592 593 594 SUBROUTINE READFTN34(NATR, IMODE4) 595C THIS SUBROUTINE READ THE ftn34 596 IMPLICIT REAL*8 (A-H,O-Z) 597 CHARACTER*80 CDUM,AGRP 598 REAL*8 DV(3,3),SOP(3,3,48),TRX(48),TRY(48),TRZ(48) 599 INTEGER C2I 600 COMMON/FTN34/ SOP,TRX,TRY,TRZ,DV,IDIM,IGROUP,NSYMOP 601 COMMON/BOHRR/ BOHR,IMODE3 602 603 include 'mode.inc' 604 605c --- DEBUG_BEGIN --- 606c print *,'in READFTN34; IMODE4=',IMODE4 607c --- DEBUG_END --- 608 609 if (IMODE4 .eq. M4_CR95) then 610 READ(34,*) CDUM 611 READ(34,*) IDIM,IGROUP 612 elseif (IMODE4 .eq. M4_CR98) then 613 READ(34,*) IDIM,IGROUP,ITYPE 614 else 615 STOP 'unknown IMODE4, must be M4_CR95 or M4_CR98' 616 endif 617 618C TAKE CARE IF CELL OF TYPE 1 IS HEXAGONAL TYPE 619 IF(IGROUP.EQ.1) THEN 620 CALL GETARG(IARG_IGRP,AGRP) 621c agrp is '8' for HEXAGONAL systems and '9' for TRIGONAL_NOT_RHOMBOHEDRAL 622 IGROUP=C2I(AGRP) 623 ENDIF 624 WRITE(12,*) 'DIM-GROUP' 625 WRITE(12,*) IDIM,IGROUP 626C VECTOR1 GOES TO FIRST COLOUMN, VETCOR2 TO SECOND & VECTOR3 TO 627C THIRD COLOUMN;;; DV(COL,ROW) 628C |11 21 31| |x1 x2 x3| 629C VEC = |12 22 32| => DV = |y1 y2 y3| 630C |13 23 33| |z1 z2 z3| 631 READ(34,*) ((DV(I,J),J=1,3),I=1,3) 632 IF(IMODE3.EQ.M3_BOHR)THEN !convert from bohrs to angstroms 633 DO I=1,3 634 DO J=1,3 635 DV(J,I) = BOHR * DV(J,I) 636 ENDDO 637 ENDDO 638 ENDIF 639 if(idim.lt.3)dv(3,3) = 1.0d0 640 if(idim.lt.2)dv(2,2) = 1.0d0 641 if(idim.lt.1)dv(1,1) = 1.0d0 642 WRITE(12,*) 'PRIMVEC' 643 644 DO I=1,3 645 WRITE(12,'(3(1x,f15.10)/3(1x,f15.10)/3(1x,f15.10))') 646 $ (DV(I,J),J=1,3) 647 ENDDO 648 READ(34,*) NSYMOP 649 DO I=1,NSYMOP 650C READ SYMM. OPERATOR 651 DO J=1,3 652 READ(34,*) (SOP(JJ,J,I),JJ=1,3) 653 ENDDO 654C READ CARTESIAN COMP. OF THE TRANSALTIUON DUE TO SYM. OPER. 655 READ(34,*) TRX(I),TRY(I),TRZ(I) 656 IF(IMODE3.EQ.M3_BOHR)THEN !convert from bohrs to angstroms 657 trx(i)=trx(i)*BOHR 658 try(i)=try(i)*BOHR 659 trz(i)=trz(i)*BOHR 660 ENDIF 661 ENDDO 662 663C READ NATR; 664 READ(34,*) NATR 665C NAT,X,Y,Z WILL ARE READ in CONVATOM subroutine !!!!!! 666 667c print *,'READFTN34_END' 668 669 RETURN 670 END 671 672 SUBROUTINE CONVATOM(IGROUP,NARG,NCSTMC,IMODE4) 673 IMPLICIT REAL*8 (A-H,O-Z) 674 include 'param.inc' 675 CHARACTER*80 FILE2 676 INTEGER NAT(NAC) 677 REAL*8 AC(3,2),BC(3,2),CC(3,2),FC(3,4),IC(3,2),RC(3,3),HC(3,3), 678 $ CSTMC(3,4), 679 * DVC(3,3),SOP(3,3,48),TRX(48),TRY(48),TRZ(48), 680 * XC(NAC,4),YC(NAC,4),ZC(NAC,4), 681 * X(NAC),Y(NAC),Z(NAC),DVEC(3,3) 682 COMMON/APOS/ AC,BC,CC,FC,IC,HC,RC,DVC,CSTMC 683 COMMON/MULTAT/ X,Y,Z,NAT,DVEC,NATR 684 COMMON/BOHRR/ BOHR,IMODE3 685 686 include 'mode.inc' 687 688c DEBUG_BEGIN 689c print *,'DEBUG: CONVATOM' 690c print *,'DEBUG: NCSTMC::', NCSTMC 691c do i=1,4 692c print *,'DEBUG: CSTMC::',(CSTMC(j,i),j=1,3) 693c enddo 694c DEBUG_END 695 696C READ NAT,X,Y,Z FROM UNIT 11 (FILE1) & PREPARE ATOMS FOR DIFFERENT 697C CENTERED CELL 698 699c print *,'IN CONVATOM' 700 701 DO I=1,NATR 702 if (narg.eq.8) then 703 READ(34,*) NAT(I),X(I),Y(I),Z(I) 704 IF(IMODE3.EQ.M3_BOHR) THEN 705 X(I) = X(I) * BOHR 706 Y(I) = Y(I) * BOHR 707 Z(I) = Z(I) * BOHR 708 ENDIF 709 endif 710 enddo 711 712 if (IMODE4 .eq. M4_CR98) then 713c generate all equivalent atoms from non-equivalent 714 call GenEquPos !GenEquPos is still bugy 715 endif 716 717c write PRIMCOOR just if we have read data from unit34 718 if (narg.eq.8) then 719 WRITE(12,*) 'PRIMCOORD' 720 WRITE(12,*) NATR,' 1' 721 endif 722 723 DO I=1,NATR 724 if(narg.eq.8) WRITE(12,'(i3,3x,3(f15.10,2x))') 725 $ NAT(I),X(I),Y(I),Z(I) 726c print *,'i;nat,x,y,z>>',i,NAT(I),X(I),Y(I),Z(I) 727C ALL DIFFERENT GROUPS HAVE (0,0,0) 728 XC(I,1)=X(I) 729 YC(I,1)=Y(I) 730 ZC(I,1)=Z(I) 731 GOTO(101,201,301,401,501,601,701,801,801,1001) IGROUP 732 733 201 CONTINUE 734C *** P->A 735 CALL GETCCOOR(DVC,AC,2,X,Y,Z,XC,YC,ZC,I,NATR) 736 GOTO 1001 737 738 301 CONTINUE 739C *** P->B 740 CALL GETCCOOR(DVC,BC,2,X,Y,Z,XC,YC,ZC,I,NATR) 741 GOTO 1001 742 743 401 CONTINUE 744C *** P->C 745 CALL GETCCOOR(DVC,CC,2,X,Y,Z,XC,YC,ZC,I,NATR) 746 GOTO 1001 747 748 501 CONTINUE 749C *** P->F 750c print *,'before DO' 751 DO IN=2,4 752 CALL GETCCOOR(DVC,FC,IN,X,Y,Z,XC,YC,ZC,I,NATR) 753 ENDDO 754c print *,'after DO' 755 GOTO 1001 756 757 601 CONTINUE 758C *** P->I 759 CALL GETCCOOR(DVC,IC,2,X,Y,Z,XC,YC,ZC,I,NATR) 760 GOTO 1001 761 762 701 CONTINUE 763C *** PR->R; primitiv rhombohedral --> rhombohedraly centered 764 DO IN=2,3 765c print *,'RC>>>>',(rc(j,in),j=1,3) 766 CALL GETCCOOR(DVC,RC,IN,X,Y,Z,XC,YC,ZC,I,NATR) 767 ENDDO 768 GOTO 1001 769 770C *** P->H 771 801 CONTINUE 772 DO IN=2,3 773 CALL GETCCOOR(DVC,HC,IN,X,Y,Z,XC,YC,ZC,I,NATR) 774 ENDDO 775C *** P->CSTMC 776 1001 CONTINUE 777 DO IN=1,NCSTMC 778 CALL GETCCOOR(DVC,CSTMC,IN,X,Y,Z,XC,YC,ZC,I,NATR) 779 ENDDO 780 101 CONTINUE 781 ENDDO 782 783C WRITE COORD. OF ATOMS WITHIN CONV. CELL IN XYZ COORD. 784 WRITE(12,*) 'CONVCOORD' 785C NCELL.......N. of ATOMS PER CONV. CELL 786 if (igroup.lt.10) then 787 NCELL=NCDET(IGROUP) 788 else 789 NCELL=NCSTMC ! ncstmc was determined by GetCenteredPoints 790 endif 791c print *,'DEBUG: NCELL=',ncell 792 WRITE(12,*) NATR,NCELL 793 DO I=1,NATR 794 DO II=1,NCELL 795 WRITE(12,'(i3,3x,3(f15.10,2x))') 796 $ NAT(I),XC(I,II),YC(I,II),ZC(I,II) 797 ENDDO 798 ENDDO 799 800c print *,'CONVATOM_END' 801 802 RETURN 803 END 804 805 806 SUBROUTINE GenEquPos 807 IMPLICIT REAL*8 (a-h,o-z) 808 include 'param.inc' 809 INTEGER NAT(NAC) 810 REAL*8 DV(3,3),SOP(3,3,48),TRX(48),TRY(48),TRZ(48), 811 * X(NAC),Y(NAC),Z(NAC),DVEC(3,3) 812 REAL*8 v(3), vr(3), dvinv(3,3), a(3), d(3) 813 REAL*8 rx(NAC), ry(NAC), rz(NAC) 814 LOGICAL not_equal, equal 815 COMMON/FTN34/ SOP,TRX,TRY,TRZ,DV,IDIM,IGROUP,NSYMOP 816 COMMON/MULTAT/ X,Y,Z,NAT,DVEC,NATR 817 818 PARAMETER (TOLMIN=1.0d-6) 819 820 call InvertSqMat123(dv, dvinv, 3) 821c print *,'DV::',((dv(i,j),i=1,3),j=1,3) 822c print *,'IDV::',((dvinv(i,j),i=1,3),j=1,3) 823 824c to fractional 825 do i=1,natr 826 vr(1)=x(i) 827 vr(2)=y(i) 828 vr(3)=z(i) 829 call ZeroMat(a,3,1) 830 call MatMult(dvinv, 3, 3, vr, 1, a) 831 rx(i)=a(1) 832 ry(i)=a(2) 833 rz(i)=a(3) 834 enddo 835 836 n_non_a=natr 837 do i=1,n_non_a 838 na = nat(i) 839 v(1) = x(i) 840 v(2) = y(i) 841 v(3) = z(i) 842 do j=1,nsymop 843 vr(1) = trx(j) 844 vr(2) = try(j) 845 vr(3) = trz(j) 846 call MatMult(sop(1,1,j),3,3,v,1,vr) 847c print *,i,j,(vr(k),k=1,3) 848c to fractional 849 call ZeroMat(a,3,1) 850 call MatMult(dvinv, 3, 3, vr, 1, a) 851c to [-0.5,0.5] 852 do i3=1,idim 853 a(i3)=a(i3) - dnint(a(i3)) 854 enddo 855c to Cartesian 856 call ZeroMat(vr,3,1) 857 call MatMult(dv, 3, 3, a, 1, vr) 858 859c check if we have new atom 860 not_equal=.true. 861 k=1 862 do while (k .le. natr .and. not_equal) 863 if ( nat(k) .ne. na ) goto 1 864 d(1)=rx(k)-a(1) 865 d(2)=ry(k)-a(2) 866 d(3)=rz(k)-a(3) 867 do i3=1,idim 868 d(i3)=d(i3) - dnint(d(i3)) 869 enddo 870 if ( abs(d(1)) + abs(d(2)) + abs(d(3)) .lt. TOLMIN ) 871 $ not_equal=.false. 872 1 continue 873 k=k+1 874 enddo 875c print *,'vr_point:',i,na,vr(1),vr(2),vr(3),.not.not_equal 876 if (not_equal) then 877c print *,'DIFF:', abs(d(1)) + abs(d(2)) + abs(d(3)) 878 natr=natr+1 879 nat(natr) = na 880 x(natr) = vr(1) 881 y(natr) = vr(2) 882 z(natr) = vr(3) 883 rx(natr)= a(1) 884 ry(natr)= a(2) 885 rz(natr)= a(3) 886c debug--debug 887c print *,'new_position: ', a(1), a(2), a(3) 888 endif 889 enddo 890 enddo 891 892 return 893 END 894 895 896 INTEGER FUNCTION NCDET(IGROUP) 897C DETERMINE HOW MENY ATOMS WILL BE IN CONV. CELL 898 IF(IGROUP.EQ.1 .OR. IGROUP.EQ.9) NCDET=1 899 IF((IGROUP.GT.1.AND.IGROUP.LT.5).OR.(IGROUP.EQ.6)) NCDET=2 900 IF(IGROUP.EQ.5) NCDET=4 901 IF(IGROUP.EQ.7.OR.IGROUP.EQ.8) NCDET=3 902 RETURN 903 END 904 905 906 SUBROUTINE GETCCOOR(A33,B33,BROW,X,Y,Z,XC,YC,ZC,NA,NATR) 907 include 'param.inc' 908 REAL*8 B33(3,*),A33(3,3),RA(3,3),COOR(3),X(NATR),Y(NATR),Z(NATR), 909 * XC(NAC,4),YC(NAC,4),ZC(NAC,4),RX(3) 910 INTEGER BROW 911 912C ******** 913C ATTENTION: HERE WE HAVE MULT. OFF "SAME-LAYING" ELEMENTS, 914C BECAUSE B33(*,BROW) IS ROW VECTOR INSTEAD OF COLOUMN !!!! 915C ******** 916c print *,'DVC::',((a33(i,j),j=1,3),i=1,3) 917c print *,'FC::',(b33(i,brow),i=1,3) 918 DO J=1,3 919 COOR(J)=0.0d0 920 DO I=1,3 921 COOR(J)=COOR(J)+A33(I,J)*B33(I,BROW) 922 ENDDO 923 ENDDO 924 XC(NA,BROW)=X(NA)+COOR(1) 925 YC(NA,BROW)=Y(NA)+COOR(2) 926 ZC(NA,BROW)=Z(NA)+COOR(3) 927 928c ***************************************************** 929c COORDINATES must be within [-1,1], check that !!! 930c ***************************************************** 931 call InvertSqMat123( A33, RA, 3 ) 932 do i=1,3 933 rx(i) = xc(na,brow)*ra(1,i) + yc(na,brow)*ra(2,i) + 934 $ zc(na,brow)*ra(3,i) 935 if ( rx(i) .gt. 1.0d0 ) then 936 rx(i) = rx(i) - 1.0d0 937 else if ( rx(i) .lt. -1.0d0 ) then 938 rx(i) = rx(i) + 1.0d0 939 endif 940 enddo 941c print *, 'DEBUG: fractC #',na,'::',rx(1),rx(2),rx(3) 942 xc(na,brow) = rx(1)*a33(1,1) + rx(2)*a33(2,1) + rx(3)*a33(3,1) 943 yc(na,brow) = rx(1)*a33(1,2) + rx(2)*a33(2,2) + rx(3)*a33(3,2) 944 zc(na,brow) = rx(1)*a33(1,3) + rx(2)*a33(2,3) + rx(3)*a33(3,3) 945 946 RETURN 947 END 948 949 950C SUBROUTINE READARG(NATOM,F1,F2,NX,NY,NZ) 951C COMMON/DIR/ NXDIR,NYDIR,NZDIR 952C CHARACTER*80 ANX,ANY,ANZ,FILE1,FILE2 953C INTEGER C2I,F1,F2 954C this command read: NXDIR, NYDIR, NZDIR arguments 955C THIS SUBROUTINE READ THE REST OF COMMAND LINE ARGUMENTS 956C & ESTIMATE THE NUMBER OF ATOMS THAT WILL BE PRODUCED 957C WITH MULTIPLICATION (COUNT JUST ONE ATOM PER CELL) 958C 959C CALL GETARG(F1,FILE1) 960C OPEN(UNIT=11,FILE=FILE1,STATUS='OLD') 961C CALL GETARG(F2,FILE2) 962C OPEN(UNIT=12,FILE=FILE2,STATUS='UNKNOWN') 963C CALL GETARG(NX,ANX) 964C CALL GETARG(NY,ANY) 965C CALL GETARG(NZ,ANZ) 966C NXDIR=C2I(ANX) 967C NYDIR=C2I(ANY) 968C NZDIR=C2I(ANZ) 969c print *,'anx,any,anz',anx,any,anz 970C NUMBER OF ATOMS/CELLS (ONE PER CELL) 971C NATOM=(NXDIR+2)*(NYDIR+2)*(NZDIR+1) 972c print *,'nx,ny,nz>',nxdir,nydir,nzdir 973C RETURN 974C END 975 976 977 978C THIS SUBROUTINE IS USED WHEN WE ARE DEALING WITH MOLECULE: 979C IT WRITE A COORDINATES OF MOLECULE/CLUSTER 980 SUBROUTINE MOLECULE 981 include 'param.inc' 982 INTEGER NAT(NAC) 983 REAL*8 X(NAC),Y(NAC),Z(NAC),DVEC(3,3) 984 COMMON/MULTAT/ X,Y,Z,NAT,DVEC,NATR 985 WRITE(12,*) 'ATOMS' 986 DO I=1,NATR 987 WRITE(12,21) NAT(I),X(I),Y(I),Z(I) 988 ENDDO 989 21 FORMAT(I5,3F16.10) 990 RETURN 991 END 992 993 LOGICAL Function IsReverse(CSTMC,RCR) 994 REAL*8 CSTMC(3,4), RC(3,3), RCR(3,3), TOLMIN 995 LOGICAL isequal 996 TOLMIN=1.0d-6 997c R(x/y/z,#point) 998 do ir=1,3 999 isequal=.FALSE. 1000c can we find RCR(*,ir) point among CSTMC(*,#) points 1001 do ic=1,3 1002 do j=1,3 1003 if ( abs(CSTMC(j,ic)-RCR(j,ir)) .lt. TOLMIN ) 1004 $ isequal=.TRUE. 1005 enddo 1006 enddo 1007 if (.not.isequal) then 1008c RCR(*,ir) point wasn't find -> OBVERSE SETTINGS 1009 IsReverse=.FALSE. 1010 return 1011 endif 1012 enddo 1013 IsReverse=.TRUE. 1014 return 1015 End 1016