c$DEBUG C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C PMIFST - Calculation of P_rincipal M_oments of I_nertia F_rom C molecular ST_ructure C C This program takes CART definitions or Cartesians as input. It first C calculates moments of inertia and then graphically displays the C structure in several ways. The structure can be rotated and structural C parameters can be checked. The program also offers a route for C conversion from Cartesian to connectivity information and high quality c PostScript graphics output via GLE. C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C NOTE: This is an IBM-PC bound program and it is tied to the graphics c package incorporated with Microsoft Fortran C C - molecule can have up to NATMAX atoms as declared in C relevant PARAMETER statements - currently 500 atoms C C ver. 7.III.2001 ----- Zbigniew KISIEL ----- C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C | http://info.ifpan.edu.pl/~kisiel/prospe.htm | C_________________________/-------------------------------------------------- C C Modification history: C c ?.?.76: First version written for IBM-360 positioned atoms as points on c 132 column lineprinter output, which could be joined by manually c drawn lines! Atom definition scheme lifted from the CART c programme c 76-95: Intermediate versions for various machines, operating systems and c graphics, including PDP-11 and ZX-Spectrum c 8.6.95: Addition of GLE copy of molecule in the displayed orientation c 17.11.95: Addition of atom labelling to GLE output c 10.10.97: Option of preserving input Cartesians for display c 12.03.98: Incremental mods. to GLE output, graphics, input menu + output c 31.03.99: 3D masking of bonds and atom circles in GLE output c 7.09.99: Graphics compatibility taken up to the level in V32 c 26.03.00: Option of reversal of axes for plotting c 16.04.00: Debugging+mods to GLE output c 30.08.00: Sorting out of axes in GLE output c 7.03.01: More GLE C C---------------------------------------------------------------------------- C C CART definition scheme: C C Atom NO=1 defines the origin, the line NO=1 to NO=2 defines the X axis C and atoms NO=1, NO=2, NO=3 define the XY plane. C C For a right handed coordinate system dihedral angles are defined C by viewing the configuration from the NC direction along NB-NA C axis or the NA-NB axis depending on whether atom C is bonded to C atom B or A. The angle is defined positive for a clockwise rota- C tion of atom C into the plane N-A-B. C............................................................................. C C 12.C Atomic masses, isotopic abundances and spins: C C 1.H 1.00782519 99.8850 1/2 C 2.H 2.0141022 0.01492 1 C 3.H 3.016050 1/2 C 10.B 10.0129388 19.78 3 C 11.B 11.0093053 80.22 3/2 C 12.C 12. 98.893 0 C 13.C 13.0033544 1.107 1/2 C 14.C 14.003241 0 C 14.N 14.0030744 99.6337 1 C 15.N 15.0001077 0.3663 1/2 C 16.O 15.9949150 99.759 0 C 17.O 16.999133 0.0374 5/2 C 18.O 17.999160 0.2039 0 C 19.F 18.9984046 100 1/2 C 20.Ne 19.9924405 90.92 (0) C 22.Ne 21.9913847 8.82 (0) C 28.Si 27.976929 92.21 (0) C 29.Si 28.976496 4.70 1/2 C 30.Si 29.973763 3.09 (0) C 31.P 30.973765 100 1/2 C 32.S 31.9720737 95 0 C 33.S 32.971462 0.760 3/2 C 34.S 33.967867 4.22 0 C 35.Cl 34.968851 75.529 3/2 C 37.Cl 36.965899 24.471 3/2 C 36.Ar 35.9675445 0.337 0 C 40.Ar 39.9623842 99.6 0 C 79.Br 78.918329 50.537 3/2 C 81.Br 80.916292 49.463 3/2 C 127.I 126.904470 100 5/2 C C............................................................................. C C When molecule is displayed under the rotation option, various actions C can be selected by the following keys: C C E,X - rotation about the horizontal axis C S,D - rotation about the vertical axis C W,R - rotation about the axis perpendicular to the plane of the screen C (the six keys above are case sensitive - lower case gives slow and C upper case fast rotation) C N - toggle display numbers of atoms C A - toggle display of principal axes c G - generation of PMIFST.GLE file displaying current orientation C P - parameter checks: display is switched to a mode allowing checks of C values of bond lengths and bond angles C In this mode: C - B selects calculation of bond length to be followed by two atom C numbers C - A selects calculation of bond angle to be followed by three atom C numbers, a,b or c can be input instead of the third number for C calculation of angle between bond and a principal axis C - D selects calculation of dihedral angle to be followed by C four atom numbers C - E exits back to previous screen C followed by Y - termination of graphics display C_____________________________________________________________________________ C C C...Initialization commands for graphics (FGRAPH.FI should only be included C into the .MAIN block, FGRAPH.FD should be declared in every program C block calling graphics routines). The three structured C variables contain coordinates: C curpos.row and curpos.col - cursor coordinates (INTEGER*2) C ixy.xcoord and ixy.ycoord - pixel coordinates (INTEGER*2) C wxy.wx and wxy.wy - window coordinates (REAL*8) C c Uncomment the statements below (and in similar headers in routines c as necessary) MSF5 = MS Fortran 5.0 c PS1 and PS4 = MS Powerstation Fortrans 1.0 and 4.0 c For MSF5 split the source before routine NFIRST c c INCLUDE 'FGRAPH.FI' MSF5+PS1 INCLUDE 'FLIB.FI' PS1 INCLUDE 'FGRAPH.FD' MSF5+PS1 c USE MSFLIB PS4 c RECORD /rccoord/curpos RECORD /xycoord/ixy RECORD /wxycoord/wxy logical*2 true C C...The code in the two INCLUDE files clashes with IMPLICIT statements - so C when using MSF5.0 graphics the recomendation is to type everything C explicitly (check with -4Yg compilation switch) C PARAMETER (NATMAX=500,NIBOND=4*NATMAX,NBONMX=6,NISOT=29, * true=.true.) COMMON /BLK0/INITIA,NATOMS,ATOMS,VBA,TRANS,PRC,RJA,VCA * /BLK1/NATMS,CORD,CORPLT,AXEND,ASPECT/BLK2/IBOND,NBONDS,W,radius COMMON /PARBLK/RANGEA,RANGEB,AMIN,BMIN COMMON /limits/maxx,maxy,linofs,curpos,ixy,wxy, * mymode,myrows,mycols,blue,red COMMON /lines/toplin,botlin,emplin,COMENT common /atnams/topicm,dtnam REAL*8 CORD(NATMAX,3),W(NATMAX),VBA(3),VCA(3),TRANS(3,3),RJA(3), * PRC(3),CORDX(NATMAX),CORDY(NATMAX),CORDZ(NATMAX),A(6), * CORPLT(NATMAX,3),corinp(natmax,3), * EIGVEC(9),ATOMS(NATMAX,3),EIGINV(3,3),AXEND(6,3) REAL*8 CON,ASYMK,BPLUSC,CS,FMX,FMY,FMZ,RANGEA,RANGEB,RANGEC, * RCA,RCB,RCC,SMALLA,SMALLB,SMALLC,SS,WW,ASPECT, * AMIN,BMIN,amax,bmax,TOPICM(NISOT),radius(natmax) INTEGER*2 maxx,maxy,linofs,mymode,myrows,mycols INTEGER*4 dummy4,blue,red INTEGER*2 IBOND(NIBOND,2),NATOMS(NATMAX,4),INKEY,N, * ICONEC(NATMAX,NBONMX),ICONN(NATMAX,2) CHARACTER*79 toplin,botlin,emplin,COMENT,LINE CHARACTER FILNAM*25,WCHAR*15,CHUNIT*8 CHARACTER*4 PMINAM(NISOT) CHARACTER*2 DTNAM(NISOT) C DATA PMINAM * /' 1H', ' 2H', ' 3H', ' 12C', * ' 13C', '29Si', ' 14N', ' 15N', * ' 16O', ' 17O', ' 18O', ' 19F', * ' 32S', ' 34S', '35Cl', '37Cl', * '79Br', '81Br', '127I', ' 6Li', * ' 7Li', ' 9Be', ' 10B', ' 11B', * '28Si', ' 31P', '36Ar', '40Ar', * ' *'/ DATA TOPICM * /1.00782519D0, 2.0141022D0, 3.016050D0, 12.0000000D0, * 13.0033544D0, 28.976496D0, 14.0030744D0, 15.0001077D0, * 15.9949150D0, 16.999133D0, 17.999160D0, 18.9984046D0, * 31.9720737D0, 33.967867D0, 34.968851D0, 36.965899D0, * 78.918329D0, 80.916292D0, 126.904470D0, 6.01513D0, * 7.01601D0, 9.01219D0, 10.01294D0, 11.009305D0, * 27.97693D0, 30.973765D0, 35.9675445D0, 39.9623842D0, * 1.0D-10/ DATA DTNAM * /' H', ' H', ' H', ' C', * ' C', 'Si', ' N', ' N', * ' O', ' O', ' O', ' F', * ' S', ' S', 'Cl', 'Cl', * 'Br', 'Br', ' I', 'Li', * 'Li', 'Be', ' B', ' B', * 'Si', ' P', 'Ar', 'Ar', * ' *'/ C WRITE(emplin,'(79(1H ))') CON=0.174532925D-1 BOHR=0.5291771D0 ITSAVE=1 C C Input of data C call graphicsmode() if($BLUE.eq.#2a0000)then dummy=setvideomode($DEFAULTMODE) call clearscreen($GCLEARSCREEN) endif WRITE(*,3000) 3000 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | PMIFST - Principal Moments of Inertia From STructural ', * 'parameters',T79,'|'/ * ' |',76(1H_),'|'/' version 7.III.2001',T64,'Zbigniew KISIEL'/) c dummy=settextcolor(11) if($BLUE.ne.#2a0000)then write(*,'(1x,a/)') * 'This program should be run in a maximised window' else write(*,'(1x,a/)') * 'This program should be run in a full screen MS-DOS window' endif dummy=settextcolor(15) C 6018 WRITE(*,6021) 6021 FORMAT(1X/' TYPE OF INPUT:'// 1 12x,' 0 = CART definitions'/ 1 12x,' 1 = Cartesians in Angstroms'/ 1 12x,'-1 = Cartesians in bohr'/ 1 12x,' 2 = CART data to be created ... ',$) READ(*,3003,ERR=6018)ICFLAG if(icflag.lt.-1.or.iflag.gt.2)goto 6018 3003 FORMAT(I5) if(icflag.eq.2)goto 3004 IF(iabs(ICFLAG).EQ.1)WRITE(*,6030) 6030 FORMAT(1X/1X,78(1H-)/' The data file should be in the' 1 ' following format:'//' line.1: comment'/ 1 ' line.2: NATOMS'/ 1 ' line.3+: N, X, Y, Z, MASS (M in amu)' 1 //' N is a four character alphanumeric atom descriptor, which is' 1 /' ignored by the program and the atoms are assigned' 1 /' consecutive numbers on input.'/1x,78(1H-)/) C 6001 WRITE(*, 3040) 3040 FORMAT(1X/' NAME OF DATA FILE : ',$) READ(*,'(A)',ERR=6001)FILNAM OPEN(2,FILE=FILNAM,STATUS='OLD',ERR=6001) READ(2,3007,err=7003)COMENT READ(2,3003,err=7003)NATMS IF(iabs(ICFLAG).EQ.1)GOTO 6019 C C...CART input from file: the card is to be prepared according to FORMAT no. C 3201; the first four integers are the defining atoms, followd by defining C bond length, bond angle, dihedral angle and mass of the atom. It is C possible to use a mnemonic for mass of the atom (as contained in PMINAM), C this should be placed in columns 60:63 of the card - just before the C position of the decimal point in a numerical mass. C DO 3206 N=1,NATMS READ(2,'(A)')LINE BACKSPACE(2) DO 4500 NN=1,NISOT IF(LINE(60:63).EQ.PMINAM(NN))THEN W(N)=TOPICM(NN) READ(2,3201,err=7004)(NATOMS(N,NNN),NNN=1,4), * (ATOMS(N,NNN),NNN=1,3) GOTO 3206 ENDIF 4500 CONTINUE READ(2,3201,err=7004)(NATOMS(N,NN),NN=1,4),(ATOMS(N,NN),NN=1,3), 1 W(N) 3206 CONTINUE 3201 FORMAT(4I4,3F13.6,F16.7) CLOSE(2) GO TO 3005 7003 write(*,'('' Error in comment or number of atoms''/)') stop 7004 write(*,'('' Error on reading atom'',i3,a/)')n,char(7) stop C C...Alternative input of Cartesian coordinates (in Angstroms or bohr) and C atomic masses (in amu). The first four characters of the line are ignored C and free format is allowed. C 6019 DO 6020 N=1,NATMS READ(2,'(A)')LINE LINE=LINE(5:) CALL NFIRST(LINE,2,CORDX(N),IDUMMY,IERROR) CALL NFIRST(LINE,2,CORDY(N),IDUMMY,IERROR) CALL NFIRST(LINE,2,CORDZ(N),IDUMMY,IERROR) CALL NFIRST(LINE,2, W(N),IDUMMY,IERROR) if(icflag.eq.-1)then cordx(n)=cordx(n)*0.5291771d0 cordy(n)=cordy(n)*0.5291771d0 cordz(n)=cordz(n)*0.5291771d0 endif corinp(n,1)=cordx(n) corinp(n,2)=cordy(n) corinp(n,3)=cordz(n) 6020 CONTINUE CLOSE(2) c c...option of avoiding inertial coordinates so that angular relationships c between structure and input axes can be followed c if(iabs(icflag).eq.1)then 7010 write(*,'(1x/'' Do you want to use input rather than '', * ''principal coordinates for graphics ? '',$)') read(*,3003,err=7010)ikeepc if(ikeepc.ne.1.and.ikeepc.ne.0)goto 7010 if(ikeepc.eq.1)then do 7011 n=1,natms cord(n,1)=cordx(n) cord(n,2)=cordy(n) cord(n,3)=cordz(n) 7011 continue goto 4010 endif endif GOTO 3005 C C...interactive CART input from keyboard C 3004 WRITE(*, 3006) 3006 FORMAT(1X/1X,'COMMENT DESCRIBING THE DATA :') READ(*,3007,ERR=3004)COMENT 3007 FORMAT(A) 6003 WRITE(*, 3008) 3008 FORMAT(1X/1X,'NUMBER OF ATOMS IN THE MOLECULE = ',$) READ(*,3003,ERR=6003)NATMS c 4530 if(iabs(icflag).ne.1)then WRITE(*,4520) 4520 FORMAT(1X/' OPTIONS: NO = 0 listing of data'/ * ' N definition of atom N'/ * ' -N change mass of atom N'/ * ' 555 backup of data') else WRITE(*,4540) 4540 FORMAT(1X/' OPTIONS: NO = 0 listing of data'/ * ' N change mass of atom N') endif c 3009 if(iabs(icflag).ne.1)then WRITE(*, 3010) 3010 FORMAT(1X/1X,'ATOM TO BE DEFINED (555=backup) NO = ',$) else write(*,4542) 4542 FORMAT(1X/1X,'ATOM whose mass is to be changed = ',$) endif READ(*,3003,ERR=3009)NO IF(NO.EQ.555)then iflag=555 GOTO 521 ENDIF IF(NO.EQ.0)GO TO 3005 if(iabs(icflag).eq.1)no=-no IF(NO.LT.0)THEN NO=-NO IF(NO.GT.NATMS)GOTO 530 GOTO 6010 ENDIF IF(NO.GT.NATMS)GOTO 530 NATOMS(NO,1)=NO c 6004 WRITE(*, 3011) 3011 FORMAT(1X, 'DEFINING ATOMS NA = ',$) READ(*,3003,ERR=6004)NATOMS(NO,2) IF(NATOMS(NO,2).LT.0.OR.NATOMS(NO,2).GE.NO)GOTO 530 c 6005 WRITE(*, 3012) 3012 FORMAT(1H+,34X,'NB = ',$) READ(*,3003,ERR=6005)NATOMS(NO,3) IF(NATOMS(NO,3).LT.0.OR.NATOMS(NO,3).GE.NO)GOTO 530 c 6006 WRITE(*, 3013) 3013 FORMAT(1H+,34X,'NC = ',$) READ(*,3003,ERR=6006)NATOMS(NO,4) IF(NATOMS(NO,4).LT.0.OR.NATOMS(NO,4).GE.NO)GOTO 530 c 6007 WRITE(*, 3014) 3014 FORMAT(1X/8X,'BOND LENGTH NO.NA = ',$) READ(*,3015,ERR=6007)ATOMS(NO,1) 3015 FORMAT(F15.10) c 6008 WRITE(*, 3016) 3016 FORMAT(8X, 'BOND ANGLE NO.NA.NB = ',$) READ(*,3015,ERR=6008)ATOMS(NO,2) c 6009 WRITE(*, 3017) 3017 FORMAT(8X, 'DIHEDRAL ANGLE NO.NA.NB.NC = ',$) READ(*,3015,ERR=6009)ATOMS(NO,3) c 6010 WRITE(*, 3018)NO 3018 FORMAT(1X/1X,'MASS OF ATOM ',I2,' (number or right justified'/ * 13x,'four character mnemonic) = ',$) READ(*,'(a)',ERR=4511)WCHAR do 4510 n=1,nisot IF(WCHAR(1:4).EQ.PMINAM(N))THEN W(NO)=TOPICM(N) goto 3009 ENDIF 4510 CONTINUE read(wchar,'(f15.9)',err=4511)w(no) GO TO 3009 4511 write(*,'('' Unable to recognise the mass, try again:'',a)') * char(7) goto 6010 c 530 WRITE(*,'(1X/'' **** ERROR IN INPUT DATA, TRY AGAIN :'',A1)') 1 CHAR(7) GOTO 3009 c 3005 if(iabs(icflag).ne.1)then WRITE(*, 3021)COMENT,NATMS 3021 FORMAT(1X/1x,78(1H-)/1x,A/1x,78(1H-)/ 1 1X,'NUMBER OF ATOMS = ',I2/1X/1X,' NO', 2 ' NA NB NC NO.NA NO.NA.NB NO.NA.NB.NC', 3 ' MASS'/) DO 3022 N=1,NATMS WRITE(*, 3023)(NATOMS(N,NN),NN=1,4),(ATOMS(N,NN),NN=1,3),W(N) 3022 CONTINUE 3023 FORMAT(1X,4I4,3F13.6,F16.7) else write(*,3046)coment,natms 3046 format(1X/1x,78(1H-)/1x,A/1x,78(1H-)/ * 1X,'NUMBER OF ATOMS = ',I2) write(*,3065) do 3047 n=1,natms write(*,3064)n,(corinp(n,i),i=1,3),w(n) 3047 continue endif 6011 WRITE(*, 3024) 3024 FORMAT(1X/' MODIFICATIONS (1/0) ? ',$) READ(*,3003,ERR=6011)IFLAG IF(IFLAG.EQ.1)GO TO 4530 if(iabs(icflag).eq.1)goto 6000 C C Initial Cartesian coordinates C CORD(1,1)=0.0D0 CORD(1,2)=0.0D0 CORD(1,3)=0.0D0 N1=1 CORD(2,1)=ATOMS(2,1) CORD(2,2)=0.0D0 CORD(2,3)=0.0D0 IF(ATOMS(3,2))5555,118,5555 118 CS=-0.3333333333D0 SS=0.94280907D0 GO TO 121 5555 CS=DCOS(CON*ATOMS(3,2)) SS=DSIN(CON*ATOMS(3,2)) NA=NATOMS(3,2) 121 IF(NA-N1)128,122,128 122 CORD(3,1)=CORD(NA,1)+ATOMS(3,1)*CS GO TO 129 128 CORD(3,1)=CORD(NA,1)-ATOMS(3,1)*CS 129 CORD(3,2)=ATOMS(3,1)*SS CORD(3,3)=0.0D0 INITIA=4 C CALL CONVER C DO 5842 N=1,NATMS CORDX(N)=CORD(N,1) CORDY(N)=CORD(N,2) CORDZ(N)=CORD(N,3) 5842 CONTINUE C C Transformation to centre of mass coordinates C 6000 WW=0.0D0 DO 991 I=1,NATMS WW=WW+W(I) 991 CONTINUE FMX=0.0D0 FMY=0.0D0 FMZ=0.0D0 DO 992 N=1,NATMS FMX=FMX+W(N)*CORDX(N) FMY=FMY+W(N)*CORDY(N) FMZ=FMZ+W(N)*CORDZ(N) 992 CONTINUE DO 993 N=1,NATMS CORDX(N)=CORDX(N)-FMX/WW CORDY(N)=CORDY(N)-FMY/WW CORDZ(N)=CORDZ(N)-FMZ/WW 993 CONTINUE C C Principal moments of inertia and eigenvectors C DO 994 N=1,6 A(N)=0.0D0 994 CONTINUE DO 995 N=1,NATMS A(1)=A(1)+W(N)*(CORDY(N)**2+CORDZ(N)**2) A(2)=A(2)-W(N)*CORDX(N)*CORDY(N) A(3)=A(3)+W(N)*(CORDX(N)**2+CORDZ(N)**2) A(4)=A(4)-W(N)*CORDX(N)*CORDZ(N) A(5)=A(5)-W(N)*CORDY(N)*CORDZ(N) A(6)=A(6)+W(N)*(CORDX(N)**2+CORDY(N)**2) 995 CONTINUE C CALL EIGEN(A,EIGVEC,3,0) C RCA=505379.07D0/A(6) RCB=505379.07D0/A(3) RCC=505379.07D0/A(1) ASYMK=(2.0*RCB-RCA-RCC)/(RCA-RCC) BPLUSC=RCB+RCC BMINC=RCB-RCC DEFECT=A(1)-A(3)-A(6) PA=0.5D0*(A(1)+A(3)-A(6)) PB=0.5D0*(A(1)+A(6)-A(3)) PC=0.5D0*(A(6)+A(3)-A(1)) C WRITE(*, 3030)RCA,A(6),PA, RCB,A(3),PB, RCC,A(1),PC, * ASYMK, DEFECT, BPLUSC, BMINC 3030 FORMAT(1X/' Rotational constants/MHz and moments of inertia/', * '(amu A**2):'// * ' A =',F13.5,' I.a = ',F12.6,' P.a = ',F12.6/ * ' B =',F13.5,' I.b = ',F12.6,' P.b = ',F12.6/ * ' C =',F13.5,' I.c = ',F12.6,' P.c = ',F12.6// * ' kappa =',F13.5,' Ic-Ia-Ib = ',F12.6// * ' B+C =',F13.5/ * ' B-C =',F13.5) C C Principal coordinates C EIGINV(1,1)=EIGVEC(5)*EIGVEC(3)-EIGVEC(6)*EIGVEC(2) EIGINV(1,2)=EIGVEC(6)*EIGVEC(1)-EIGVEC(4)*EIGVEC(3) EIGINV(1,3)=EIGVEC(4)*EIGVEC(2)-EIGVEC(5)*EIGVEC(1) EIGINV(2,1)=EIGVEC(2)*EIGVEC(9)-EIGVEC(8)*EIGVEC(3) EIGINV(2,2)=EIGVEC(7)*EIGVEC(3)-EIGVEC(9)*EIGVEC(1) EIGINV(2,3)=EIGVEC(8)*EIGVEC(1)-EIGVEC(7)*EIGVEC(2) EIGINV(3,1)=EIGVEC(8)*EIGVEC(6)-EIGVEC(5)*EIGVEC(9) EIGINV(3,2)=EIGVEC(9)*EIGVEC(4)-EIGVEC(7)*EIGVEC(6) EIGINV(3,3)=EIGVEC(7)*EIGVEC(5)-EIGVEC(8)*EIGVEC(4) WW=EIGVEC(7)*EIGINV(1,1)+EIGVEC(4)*EIGINV(2,1) 1 +EIGVEC(1)*EIGINV(3,1) DO 3060 N=1,3 DO 3060 NN=1,3 EIGINV(N,NN)=EIGINV(N,NN)/WW 3060 CONTINUE DO 997 N=1,NATMS CORD(N,1)=CORDX(N)*EIGINV(1,1)+CORDY(N)*EIGINV(1,2)+ 1 CORDZ(N)*EIGINV(1,3) CORD(N,2)=CORDX(N)*EIGINV(2,1)+CORDY(N)*EIGINV(2,2)+ 1 CORDZ(N)*EIGINV(2,3) CORD(N,3)=CORDX(N)*EIGINV(3,1)+CORDY(N)*EIGINV(3,2)+ 1 CORDZ(N)*EIGINV(3,3) 997 CONTINUE C C Hardcopy of results C 6012 WRITE(*, 3061) 3061 FORMAT(1X/' Principal coordinates and hardcopy of results ?'/ * ' (-1 for coordinates in bohr) ... ' * ,$) READ(*,3003,ERR=6012)IFLAG IF(IABS(IFLAG).NE.1)GO TO 4010 CHUNIT='Angstrom' IF(IFLAG.EQ.-1)CHUNIT='Bohr' C OPEN(3,FILE='PMIFST.RES',STATUS='UNKNOWN',ERR=6016) GOTO 6017 6016 WRITE(*,'(1X/'' **** ERROR on creation of output file - '', 1 '' try a different file name: '',$)') READ(*,'(A)',ERR=6016)FILNAM OPEN(3,FILE=FILNAM,STATUS='UNKNOWN',ERR=6016) C 6017 WRITE(3,2002)COMENT 2002 FORMAT(1X,78(1H-)/1x,A/1X,78(1H-)//' Input data:'/) c if(iabs(icflag).eq.1)then write(3,3065) do 3066 n=1,natms write(3,3064)n,(corinp(n,i),i=1,3),w(n) 3066 continue else WRITE(3,3067) 3067 FORMAT(1X,' NO', 1 ' NA NB NC NO.NA NO.NA.NB NO.NA.NB.NC', 2 ' MASS'/) DO 3068 N=1,NATMS WRITE(3, 3023)(NATOMS(N,NN),NN=1,4),(ATOMS(N,NN),NN=1,3),W(N) 3068 CONTINUE endif 3065 FORMAT(1X,'ATOM NO.',7X,'x',14X,'y',14X,'z',13X,'MASS'/) c DO 2000 N=1,2 NOUT=0 IF(N.EQ.2)NOUT=3 c if(n.eq.2)then write(nout,'(1x/'' Rotation matrix:''/10x,''a'',12x,''b'', * 12x,''c'')') write(nout,'(1x,a1,3f13.8)')'x',eigvec(7),eigvec(8),eigvec(9) write(nout,'(1x,a1,3f13.8)')'y',eigvec(4),eigvec(5),eigvec(6) write(nout,'(1x,a1,3f13.8)')'z',eigvec(1),eigvec(2),eigvec(3) endif c WRITE(NOUT,3063)chunit 3063 FORMAT(1X/' Principal coordinates /',a,':'// * 1X,'ATOM NO.',7X,'A',14X,'B',14X,'C',13X,'MASS'/) DO 2001 I=1,NATMS IF(IFLAG.EQ.1)THEN WRITE(NOUT,3064)I,CORD(I,1),CORD(I,2),CORD(I,3),W(I) ELSE WRITE(NOUT,3064)I,CORD(I,1)/BOHR,CORD(I,2)/BOHR, * CORD(I,3)/BOHR,W(I) ENDIF 2001 CONTINUE 3064 FORMAT(2X,I2,2X,3F15.6,F17.7) 2000 CONTINUE WRITE(3, 3030)RCA,A(6),PA, RCB,A(3),PB, RCC,A(1),PC, * ASYMK, DEFECT, BPLUSC, BMINC WRITE(*,'(/'' Results have been written to PMIFST.RES'',A1)') 1 CHAR(7) write(3,3069) 3069 FORMAT(1x/1X,78(1H-)) CLOSE(3) C C...Location of bonds, determination of scaling parameters and of points C defining principal axes C 4010 CALL BONLOC C DO 3500 N=1,NATMS CORPLT(N,1)=CORD(N,1) CORPLT(N,2)=CORD(N,2) CORPLT(N,3)=CORD(N,3) 3500 CONTINUE C DO 510 N=1,6 DO 510 NN=1,3 AXEND(N,NN)=0.D0 510 CONTINUE CALL SCALIN(1,1,RANGEA,SMALLA) AXEND(1,1)=SMALLA AXEND(2,1)=SMALLA+RANGEA CALL SCALIN(2,1,RANGEB,SMALLB) AXEND(3,2)=SMALLB AXEND(4,2)=SMALLB+RANGEB CALL SCALIN(3,1,RANGEC,SMALLC) AXEND(5,3)=SMALLC AXEND(6,3)=SMALLC+RANGEC IF(RANGEC.LT.2.D0)THEN AXEND(5,3)=-1.D0 AXEND(6,3)=1.D0 ENDIF C C C...Plots of projections on the three principal planes C 6014 if(ikeepc.ne.1)then WRITE(*, 4001)'principal' else WRITE(*, 4001)'Cartesian' endif 4001 FORMAT(1X/' Projections on the three ',a,' planes ? ',$) READ(*,3003,ERR=6014)IFLAG IF(IFLAG.NE.1)GO TO 4023 C C...initialise graphics - see comments in routine ROTATE for details call graphicsmode() if(myrows.eq.25.or.myrows.eq.30)then mult=1 else mult=2 endif ASPECT=maxy/(0.75*maxx) dummy4 = setbkcolor( BLUE ) pixcol=getcolor() $NODEBUG CALL setlinestyle(#FFFF) $DEBUG CALL moveto(0,linofs,ixy) dummy=lineto(maxx,linofs) CALL moveto(0,int2(maxy-mult*linofs),ixy) dummy=lineto(maxx,maxy-mult*linofs) WRITE(toplin,'(a)')coment WRITE(botlin,'(2A)')'Projections on the three principal planes', 1 ' Press ENTER to exit' if(ikeepc.eq.1)botlin(26:34)='Cartesian' CALL settextposition(1,1,curpos) CALL outtext(toplin) CALL settextposition(myrows,1,curpos) CALL outtext(botlin) c IH=maxx/2 IV=maxy/2 C C...A-B plane CALL SCALIN(1,1,RANGEA,SMALLA) CALL SCALIN(2,1,RANGEB,SMALLB) AMEAN=SMALLA+0.5D0*RANGEA BMEAN=SMALLB+0.5D0*RANGEB RANGEA=1.05D0*RANGEA RANGEB=1.05D0*RANGEB PIXELR=REAL(IV)/REAL(IH) IF(RANGEB*ASPECT/RANGEA.LE.PIXELR)THEN RANGEB=RANGEA*PIXELR/ASPECT ELSE RANGEA=RANGEB*ASPECT/PIXELR ENDIF AMIN=AMEAN-RANGEA*0.5D0 AMAX=AMEAN+RANGEA*0.5D0 BMIN=BMEAN-0.5D0*RANGEB BMAX=BMEAN+0.5D0*RANGEB CALL setviewport(0,linofs+1,maxx-IH,linofs+1+IV) dummy=setwindow(TRUE,amin,bmin,amax,bmax) CALL view(1,2,1,1,1,0,ikeepc) C C...B-C plane CALL SCALIN(2,1,RANGEA,SMALLA) CALL SCALIN(3,1,RANGEB,SMALLB) AMEAN=SMALLA+0.5D0*RANGEA BMEAN=SMALLB+0.5D0*RANGEB RANGEA=1.05D0*RANGEA RANGEB=1.05D0*RANGEB IF(RANGEB*ASPECT/RANGEA.LE.PIXELR)THEN RANGEB=RANGEA*PIXELR/ASPECT ELSE RANGEA=RANGEB*ASPECT/PIXELR ENDIF AMIN=AMEAN-RANGEA*0.5D0 AMAX=AMEAN+RANGEA*0.5D0 BMIN=BMEAN-0.5D0*RANGEB BMAX=BMEAN+0.5D0*RANGEB CALL setviewport(maxx-IH,linofs+1,maxx,linofs+1+IV) dummy=setwindow(TRUE,amin,bmin,amax,bmax) CALL view(2,3,1,1,1,0,ikeepc) C C...A-C plane CALL SCALIN(1,1,RANGEA,SMALLA) CALL SCALIN(3,1,RANGEB,SMALLB) AMEAN=SMALLA+0.5D0*RANGEA BMEAN=SMALLB+0.5D0*RANGEB RANGEA=1.05D0*RANGEA RANGEB=1.05D0*RANGEB IF(RANGEB*ASPECT/RANGEA.LE.PIXELR)THEN RANGEB=RANGEA*PIXELR/ASPECT ELSE RANGEA=RANGEB*ASPECT/PIXELR ENDIF AMIN=AMEAN-RANGEA*0.5D0 AMAX=AMEAN+RANGEA*0.5D0 BMIN=BMEAN-0.5D0*RANGEB BMAX=BMEAN+0.5D0*RANGEB IHH=maxx/2-IH/2 CALL setviewport(IHH,maxy-linofs-1-IV,IHH+IH,maxy-mult*linofs-1) dummy=setwindow(TRUE,amin,bmin,amax,bmax) CALL view(1,3,1,1,1,0,ikeepc) c dummy=settextcolor(11) CALL settextposition(3,1,curpos) if(ikeepc.eq.0)then CALL outtext('A-B') else CALL outtext('X-Y') endif CALL settextposition(3,int2(mycols/2),curpos) if(ikeepc.eq.0)then CALL outtext('B-C') else CALL outtext('Y-Z') endif CALL settextposition(int2(myrows-2),1,curpos) if(ikeepc.eq.0)then CALL outtext('A-C') else CALL outtext('X-Z') endif dummy=settextcolor(15) C 500 IK=INKEY(N) IF(IK.NE.13)GOTO 500 dummy=setvideomode($DEFAULTMODE) call clearscreen($GCLEARSCREEN) C C...Rotations of the molecule and structure checks C 4023 WRITE(*,'(/'' Rotations and structure checks ? '',$)') READ(*,3003,ERR=4023)IFLAG IF(IFLAG.NE.1)GOTO 3062 CALL ROTATE(ikeepc) C C...final decisions for Cartesian input C 3062 IF(iabs(ICFLAG).EQ.1)THEN ikeepc=0 350 WRITE(*,351) 351 FORMAT(1X/ * ' OPTIONS: 0 = EXIT'/ * ' 1 = go to beginning'/ * ' 2 = produce connectivity output ... ',$) READ(*,3003,ERR=350)IFLAG IF(IFLAG.LT.0.OR.IFLAG.GT.2)GOTO 350 IF(IFLAG.EQ.0)GOTO 3043 IF(IFLAG.EQ.1)GOTO 4530 IF(IFLAG.EQ.2)GOTO 521 ENDIF C C...final decisions for CART input C WRITE(*, 3044) 3044 FORMAT(1X/' MODIFICATIONS ? ',$) READ(*,3003,ERR=3062)IFLAG IF(IFLAG.EQ.1)GO TO 4530 C C Hardcopy of data C ITSAVE=0 6015 WRITE(*, 3045) 3045 FORMAT(1X/' DATA OUTPUT 0 = none'/ * ' 1 = CART standard'/ * ' 2 = DTMM standard ..... ',$) READ(*,3003,ERR=6015)IFLAG IF(IFLAG.LT.0.OR.IFLAG.GT.2)GOTO 6015 C 521 IF(IFLAG.EQ.1.OR.IFLAG.EQ.555)THEN 520 WRITE(*, 3042) 3042 FORMAT(1X/' Name of output file : ',$) READ(*,'(A)',ERR=520)FILNAM OPEN(3,FILE=FILNAM,STATUS='UNKNOWN',ERR=520) WRITE(3,3007)COMENT WRITE(3,3003)NATMS DO 3205 N=1,NATMS WRITE(3,3201)(NATOMS(N,NN),NN=1,4),(ATOMS(N,NN),NN=1,3), 1 W(N) 3205 CONTINUE CLOSE(3) IF(ITSAVE.EQ.1)GOTO 4530 ENDIF C C...Produce output file in standard of DTMM .MOL files: C ICONEC is the connectivity matrix set up on the basis of bonds C located with BONLOC C ICONN(N,1) contains the number of atoms connected to atom N C ICONN(N,2) contains the lookup address for name of atom N C IF(IFLAG.EQ.2)THEN DO 541 N=1,NATMS DO 541 NN=1,NBONMX ICONEC(N,NN)=0 541 CONTINUE C C...work out the connectivity and lookup address C DO 538 N=1,NATMS DO 536 NN=1,NISOT IF(NINT(W(N)).EQ.NINT(TOPICM(NN)))THEN ICONN(N,2)=NN GOTO 537 ENDIF 536 CONTINUE ICONN(N,2)=NISOT C 537 NCONAT=0 DO 531 NN=1,NBONDS IF(IBOND(NN,1).EQ.IBOND(NN,2))GOTO 531 IF(IBOND(NN,1).EQ.N)THEN NCONAT=NCONAT+1 ICONEC(N,NCONAT)=IBOND(NN,2) ENDIF IF(IBOND(NN,2).EQ.N)THEN NCONAT=NCONAT+1 ICONEC(N,NCONAT)=IBOND(NN,1) ENDIF 531 CONTINUE ICONN(N,1)=NCONAT 538 CONTINUE C 532 WRITE(*,3042) READ(*,'(A)',ERR=532)FILNAM OPEN(3,FILE=FILNAM,STATUS='UNKNOWN',ERR=532) C WRITE(3,533)1.0,1.0,1.0 533 FORMAT(40X,F6.3,3X,F6.3,3X,F6.3) WRITE(3,534)90.0,90.0,90.0 534 FORMAT(22X,F7.3,1X,F7.3,1X,F7.3) WRITE(3,535)NATMS,COMENT 535 FORMAT(I4/1X,A) C DO 539 N=1,NATMS WRITE(3,540)N,DTNAM(ICONN(N,2)),CORD(N,1),CORD(N,2), * CORD(N,3),(ICONEC(N,NN),NN=1,ICONN(N,1)) 539 CONTINUE 540 FORMAT(I4,A2,4X,3(F10.5),1X,6(I4)) ENDIF C 3043 STOP END C C_____________________________________________________________________________ C SUBROUTINE graphicsmode() C C This routine sets the graphics mode on the PC and C determines the pixel limits on x and y axes (0,maxx), (0,maxy) C INCLUDE 'FGRAPH.FD' MSF5+PS1 c USE MSFLIB PS4 C RECORD /rccoord/curpos RECORD /xycoord/ixy RECORD /wxycoord/wxy INTEGER*2 dummy,maxx,maxy,linofs,mymode,myrows,mycols INTEGER*4 blue,red RECORD /videoconfig/myscreen COMMON /limits/maxx,maxy,linofs,curpos,ixy,wxy, * mymode,myrows,mycols,blue,red C C...Set graphics mode. C CALL getvideoconfig(myscreen) C C...this setting was an old default (EGA and 25 rows of text) c mymode=$ERESCOLOR c myrows=25 C C...this setting should have widest compatibility (640*480, 30 line VGA) c mymode=$VRES256COLOR mymode=$VRES16COLOR myrows=30 C C...this is for 800*600 SVGA, usually 16 colors is better than 256 c mymode=$SRES256COLOR c mymode=$SRES16COLOR c myrows=37 c dummy = setvideomoderows ( mymode, myrows ) myrows=dummy if(myrows.ne.25.and.myrows.ne.30.and.myrows.ne.37)then write(*,'(1x//i5,a,a//)')myrows, * ' <--- this is the number of lines in this mode', * ' and it is not yet supported' pause endif C IF( dummy .EQ. 0 ) STOP 'Error: cannot set graphics mode' C C...Determine parameters of the mode C CALL getvideoconfig( myscreen ) maxx = myscreen.numxpixels - 1 maxy = myscreen.numypixels - 1 mycols=myscreen.numtextcols linofs=nint(real(maxy+1)/real(myrows)) c c...patch for different use of SETBKCOLOR by PS4 from PS1 (and MSF5) c blue=$BLUE red=$RED if($BLUE.ne.#2a0000)then blue=1 red=4 endif c dummy4 = setbkcolor( BLUE ) call clearscreen($GCLEARSCREEN) C RETURN END C_____________________________________________________________________________ C SUBROUTINE VIEW(IX,IY,ixrev,iyrev,ICONTR,ICLEAR,ikeepc) C C Graphics plot of rectangular projection of molecular structure onto C one of the principal planes C C NOTE: the pixel origin is in the top lh. corner with Y increasing C downwards, but in the floating point window declared with the C option .TRUE. the Y-axis increases in the usual up direction. C C IX, IY - principal axes to be used for X and Y coordinates of C the plot (=1,2,3 for a,b,c resp.) c if IX is set negative then GLE output is produced c IXREV,IYREV - whether IX,IY are to be reversed (-1) or not (1) C ICONTR = 0 draw the molecule C +1 draw the molecule and principal axes (do not label) C +2 draw the molecule and principal axes (label axes) C +10 number atoms C ICLEAR = 0 do not clear display prior to drawing C 1 clear display C INCLUDE 'FGRAPH.FD' MSF5+PS1 c USE MSFLIB PS4 c RECORD /rccoord/curpos RECORD /xycoord/ixy RECORD /wxycoord/wxy C parameter (nisot=29) COMMON /TXTBUF/LINE COMMON /lines/toplin,botlin,emplin,COMENT COMMON /limits/maxx,maxy,linofs,curpos,ixy,wxy, * mymode,myrows,mycols,blue,red COMMON /BLK1/NATMS,CORD,CORPLT,AXEND,ASPECT/BLK2/IBOND,NBONDS,WT, * radius COMMON /PARBLK/RANGEX,RANGEY,XMIN,YMIN common /atnams/topicm,dtnam common /sort/zcord,izpt PARAMETER (NATMAX=500,NIBOND=4*NATMAX) REAL*8 CORD(NATMAX,3),WT(NATMAX),CORPLT(NATMAX,3),AXEND(6,3), * XMIN,YMIN,topicm(nisot),radius(natmax),zcord(nibond) REAL*8 ASPECT,RANGEX,RANGEY,x1,y1,x2,y2,DELX,DELY INTEGER*2 maxx,maxy,mymode,myrows,mycols,linofs,dummy INTEGER*4 blue,red INTEGER*2 IBOND(NIBOND,2) INTEGER*4 ISTART,JSTART,izpt(nibond) CHARACTER*79 toplin,botlin,emplin,COMENT CHARACTER LINE*80,AXNAM(6),dtnam*2(nisot),cdummy*1 integer*2 tday,tmon,tyear,thour,tmin,tsec,t100 DATA AXNAM/'A','B','C','X','Y','Z'/ C iaxis=icontr-(icontr/10)*10 inumbr=icontr/10 iz=iabs(6/(ix*iy)) c if(myrows.eq.25.or.myrows.eq.30)then mult=1 else mult=2 endif IYSPAN=maxy-(1+mult)*linofs if(ix.lt.0)then ix=-ix igle=1 axoffs=1.5 open(4,file='pmifst.gle',status='unknown') write(4,40) 40 format(' !',48x,'PMIFST version 7.III.2001'/' !',72(1h-)/' !') write(4,35)'! ',coment write(4,39) 39 format(' !'/' !',72(1h-)/' !') write(4,34)'size 29.8 21' write(4,34)'!' write(4,34)'! PLOT CUSTOMISATION, set the first four control ' * //'variables as follows:' write(4,34) * '! 0,1,0,1 = default, ball&stick, heavy atom labels only' write(4,34)'! 0,1,1,1 = ball&stick, all atoms labelled' write(4,34)'! 0,1,0,0 = ball&stick, no labels' write(4,34)'! 1,0,0,0 = long sticks only' write(4,34)'!' c 511 CALL clearscreen($GCLEARSCREEN) CALL settextposition(3,1,curpos) radmul=0.30 write(*,510)radmul 510 format(1x/' ----> Default ball radius as fraction of ', * 'covalent radius is',f5.2// * 8x,'Press ENTER to accept or type in new value: ',$) read(*,'(f12.8)',err=511)r1 if(r1.lt.0.0)goto 511 if(r1.gt.0.0)radmul=r1 c glesx=20.0/rangex 211 write(*,210)glesx 210 format(1x///' ----> Current scale is',f10.4,' cm/Angstr'// * 8x,'Press ENTER to accept or type in new value: ',$) read(*,'(f12.8)',err=211)r1 if(r1.lt.0.0)goto 211 if(r1.gt.0.0)glesx=r1 glesy=glesx c r1=0.95*glesx*radmul write(line,'(f10.2)')r1 call cleans(line,10,nst) write(4,34)'longst=0' write(4,34)'atball=1' write(4,34)'labelh=0' write(4,34)'labela=1' write(4,34)'labsiz='//line(nst:10) write(4,34)'labf$="texcmssb" ! use texcmr or texcmssb' write(4,39) c write(line,'(f10.5)')radmul call cleans(line,10,nst) write(4,34)'radmul='//line(nst:10) write(4,34)'hbondw=0.1' write(4,34)'tbondw=0.2' write(4,34)'circw=0.02' write(4,34)'cirlet=0.45*labsiz' cirlet=0.45 c addx=2. addy=1.5 xorig=-xmin*glesx+addx yorig=-ymin*glesy+addy write(4,34)'!' write(4,36)'! scaling =',glesx,' Angstr/cm' write(4,36)'! xorig =',xorig write(4,36)'! yorig =',yorig write(4,34)'!' write(4,34)'amove 0 0' write(4,34)'begin translate ',0.5*(29.8-2.*xorig), * 0.5*(21.-2.*yorig) goto 33 else igle=0 endif 36 format(1x,a,f10.5,a,f10.5) c 212 IF(ICLEAR.EQ.1)CALL clearscreen($GVIEWPORT) C C...Plot cirles at atom positions (small circles normally, large circles C for atom numbers) C DELX=RANGEX/50.D0 DELY=RANGEY/50.D0 DO 3 N=1,NATMS IF(inumbr.LT.1)THEN X1=ixrev*CORPLT(N,IX)-DELX*0.3D0 Y1=iyrev*CORPLT(N,IY)+DELY*0.3D0 X2=ixrev*CORPLT(N,IX)+DELX*0.3D0 Y2=iyrev*CORPLT(N,IY)-DELY*0.3D0 ELSE IF(WT(N).LT.0.1D0)GOTO 3 X1=ixrev*CORPLT(N,IX)-DELX Y1=iyrev*CORPLT(N,IY)+DELY X2=ixrev*CORPLT(N,IX)+DELX Y2=iyrev*CORPLT(N,IY)-DELY ENDIF dummy=ellipse_w($GFILLINTERIOR,X1,Y1,X2,Y2) 3 CONTINUE C C...Plot and label axes C IF(Iaxis.lt.1)GOTO 110 $NODEBUG CALL setlinestyle(#8888) $DEBUG 33 if(igle.eq.1)then write(4,39) write(4,34)'! PRINCIPAL AXES' write(4,39) write(4,34)'set lwidth 0.01' write(4,34)'set hei 0.7' write(4,34)'set font texcmmi' nn=0 endif dummy=setcolor(11) DO 4 N=1,5,2 if(igle.ne.1)then x1=ixrev*AXEND(N,IX) y1=iyrev*AXEND(N,IY) call moveto_w(x1,y1,wxy) x2=ixrev*AXEND(N+1,IX) y2=iyrev*AXEND(N+1,IY) dummy=lineto_w(x2,y2) else nn=nn+1 write(4,34)'!' write(4,34)'! '//axnam(nn+ikeepc*3)//'-axis' write(4,34)'!' xgle=(ixrev*axend(n,ix)-xmin)*glesx+addx ygle=(iyrev*axend(n,iy)-ymin)*glesy+addy xgle1=(ixrev*axend(n+1,ix)-xmin)*glesx+addx ygle1=(iyrev*axend(n+1,iy)-ymin)*glesy+addy c write(4,34)'amove ',xgle,ygle c write(4,34)'aline ',xgle1,ygle1 if(abs(xgle1-xgle).lt.0.1.and.abs(ygle1-ygle).lt.0.1)goto 4 if(abs(xgle1-xgle).gt.1.e-10)then grad=(ygle1-ygle)/(xgle1-xgle) else grad=1.E+10 endif xgle2=sign(1.,xgle1-xorig)*axoffs/dsqrt(1.d0+grad**2)+xgle1 ygle2=sign(1.,ygle1-yorig) * *abs(grad)*axoffs/dsqrt(1.d0+grad**2)+ygle1 write(4,34)'set lstyle 9' write(4,34)'amove ',xgle2,ygle2 xgle3=sign(1.,xgle-xorig)*axoffs/dsqrt(1.d0+grad**2)+xgle ygle3=sign(1.,ygle-yorig) * *abs(grad)*axoffs/dsqrt(1.d0+grad**2)+ygle write(4,34)'aline ',xgle3,ygle3 c NN=N if(ixrev*axend(n,ix) +iyrev*axend(n,iy).lt. * ixrev*axend(n+1,ix)+iyrev*axend(n+1,iy))NN=n+1 write(4,34)'set lstyle 0' xgle=(ixrev*axend(nn,ix)-xmin)*glesx+addx ygle=(iyrev*axend(nn,iy)-ymin)*glesy+addy xgle3=sign(1.,xgle-xorig) * *(axoffs+0.5)/dsqrt(1.d0+grad**2)+xgle ygle3=sign(1.,ygle-yorig) * *abs(grad)*(axoffs+0.5)/dsqrt(1.d0+grad**2)+ygle write(4,34)'amove ',xgle3-0.2,ygle3-0.2 write(4,35)'text ',char(ichar(axnam((n+1)/2+ikeepc*3))+32) c endif 4 CONTINUE dummy=setcolor(15) 34 format(1x,a,2f10.5) 35 format(1x,5a) $NODEBUG CALL setlinestyle(#FFFF) $DEBUG IF(Iaxis.eq.2)THEN dummy=setcolor(11) DO 5 N=2,6,2 NN=N if(ixrev*axend(n,ix) +iyrev*axend(n,iy).lt. * ixrev*axend(n-1,ix)+iyrev*axend(n-1,iy))NN=n-1 if(igle.ne.1)then ISTART=maxx*(ixrev*AXEND(NN,IX)-XMIN)/RANGEX-4 JSTART=IYSPAN-IYSPAN*(iyrev*AXEND(NN,IY)-YMIN)/RANGEY+4 WRITE(LINE,'(A)')AXNAM(N/2+ikeepc*3) CALL PLOTCH(2,2,ISTART,JSTART,-1) CALL PLOTCH(2,2,ISTART+1,JSTART,-1) endif 5 CONTINUE dummy=setcolor(15) ENDIF C C...Plot bonds, heavy atom bonds with continuous line, bonds to hydrogen C with dotted line, dummy atoms not bonded C 110 if(igle.ne.1)then DO 8 N=1,NBONDS IF(IBOND(N,1).EQ.IBOND(N,2))GOTO 8 C IF(WT(IBOND(N,1)).LE.2.1D0.OR.WT(IBOND(N,2)).LE.2.1D0)THEN $NODEBUG CALL setlinestyle(#AAAA) $DEBUG ELSE $NODEBUG CALL setlinestyle(#FFFF) $DEBUG ENDIF CALL moveto_w(ixrev*CORPLT(IBOND(N,1),IX), 1 iyrev*CORPLT(IBOND(N,1),IY),wxy) dummy=lineto_w(ixrev*CORPLT(IBOND(N,2),IX), 1 iyrev*CORPLT(IBOND(N,2),IY)) 8 CONTINUE endif c c...For GLE output plot sticks and atoms in order of decreasing depth c to obtain proper masking effects c if(igle.eq.1)then c c...set up depth array, +ve IZPT pointers = atoms, -ve = bonds c nz=0 do 205 n=1,natms nz=nz+1 izpt(nz)=n zcord(nz)=corplt(n,iz) 205 continue do 206 n=1,nbonds x1=ixrev*corplt(ibond(n,1),ix) y1=iyrev*corplt(ibond(n,1),iy) z1=corplt(ibond(n,1),iz) x2=ixrev*corplt(ibond(n,2),ix) y2=iyrev*corplt(ibond(n,2),iy) z2=corplt(ibond(n,2),iz) rgle=dsqrt( (x2-x1)**2 + (y2-y1)**2 + (z2-z1)**2 ) if(rgle.eq.0.0)rgle=1.E-10 r1=radius(ibond(n,1))*radmul r2=radius(ibond(n,2))*radmul z3=(z1+r1*(z2-z1)/rgle) z4=(z2-r2*(z2-z1)/rgle) nz=nz+1 izpt(nz)=-n zcord(nz)=z3 if(z4.lt.z3)zcord(nz)=z4 206 continue c c...sort pointers in IZPT in order of increasing values in ZCORD c call sortp(1,int4(nz)) c c...plot elements in the order in IZPT c write(4,39) write(4,34)'! BOND STICKS AND ATOM LABELS' write(4,39) write(4,34)'set font labf$' do 207 nnn=1,nz n=izpt(nnn) c c...bond N c if(n.lt.0)then n=-n write(4,34)'!' write(4,37)'! Bond ',ibond(n,1),' -',ibond(n,2) write(4,34)'!' 37 format(1x,a,i3,a,i3) write(4,34)'set lstyle 0' IF(WT(IBOND(N,1)).LE.2.1D0.OR.WT(IBOND(N,2)).LE.2.1D0)THEN write(4,34)'set lwidth hbondw' else write(4,34)'set lwidth tbondw' endif write(4,34)'if longst=1 then' xgle1=(ixrev*corplt(ibond(n,1),ix)-xmin)*glesx+addx ygle1=(iyrev*corplt(ibond(n,1),iy)-ymin)*glesy+addy xgle2=(ixrev*corplt(ibond(n,2),ix)-xmin)*glesx+addx ygle2=(iyrev*corplt(ibond(n,2),iy)-ymin)*glesy+addy write(4,34)' amove ',xgle1,ygle1 write(4,34)' aline ',xgle2,ygle2 write(4,34)'else' c x1=ixrev*corplt(ibond(n,1),ix) y1=iyrev*corplt(ibond(n,1),iy) z1=corplt(ibond(n,1),iz) x2=ixrev*corplt(ibond(n,2),ix) y2=iyrev*corplt(ibond(n,2),iy) z2=corplt(ibond(n,2),iz) rgle=dsqrt( (x2-x1)**2 + (y2-y1)**2 + (z2-z1)**2 ) if(rgle.eq.0.0)rgle=1.E-10 r1=radius(ibond(n,1))*radmul r2=radius(ibond(n,2))*radmul xgle3=(x1-xmin+r1*(x2-x1)/rgle)*glesx+addx ygle3=(y1-ymin+r1*(y2-y1)/rgle)*glesy+addy xgle4=(x2-xmin-r2*(x2-x1)/rgle)*glesx+addx ygle4=(y2-ymin-r2*(y2-y1)/rgle)*glesy+addy write(4,34)' amove ',xgle3,ygle3 write(4,34)' aline ',xgle4,ygle4 write(4,34)'end if' else c c...circle and label for atom N: atom labels and underlying clearing circle c are fixed in size, whereas circle from atomic sphere scales according to c size of the molecule. If sphere becomes smaller than annotation then c the latter is skipped c do 201 nn=1,nisot if( nint(wt(n)) .eq. nint(topicm(nn)) )goto 202 201 continue goto 207 202 xgle=(ixrev*corplt(n,ix)-xmin)*glesx+addx ygle=(iyrev*corplt(n,iy)-ymin)*glesy+addy c write(4,35)'!' write(4,37)'! Atom ',n write(4,34)'!' write(4,34)'set hei labsiz' write(4,34)'set lwidth circw' write(4,34)'if atball=1 then' write(4,34)' amove ',xgle,ygle rad=radius(n)*glesx write(4,334)' circle ',rad,'*radmul fill white' write(4,34)'end if' 334 format(a,f5.2,a) c cdummy='h' if(nint(wt(n)).gt.3)cdummy='a' write(4,34)'if label'//cdummy//'=1 then' write(4,34)' amove ',xgle,ygle write(4,34)' set color white' write(4,34)' circle cirlet fill white' 38 format(1x,a,f5.3,a) write(4,34)' set color black' if(dtnam(nn)(1:1).eq.' ')then write(4,35)' rmove -0.3*labsiz -0.3*labsiz' write(4,35)' text ',dtnam(nn)(2:2) else write(4,35)' rmove -0.6*labsiz -0.3*labsiz' write(4,35)'text ',dtnam(nn) endif write(4,34)'end if' c endif 207 continue c write(4,34)'!' write(4,34)'end translate' write(4,39) call getdat(tyear,tmon,tday) call gettim(thour,tmin,tsec,t100) write(4,'(a)')' set hei 0.7' write(4,'(a)')' amove 1 20.0' write(4,1500)thour,':',tmin,':',tsec,tday,'/',tmon,'/',tyear 1500 format(' text {\bf PMIFST}_{',i2,a,i2,a,i2,' \,', * i2,a,i2,a,i4,'}') write(4,40) close(4) write(*,'(1x///1x,'' File PMIFST.GLE has been generated, '', * ''press E N T E R to continue '',$)') read(*,'(i1)',err=215)igle 215 igle=0 c CALL clearscreen($GVIEWPORT) c CALL setviewport(0,linofs+1,maxx-1,maxy-mult*linofs-1) c CALL moveto(0,0,ixy) c dummy=lineto(0,maxy-(1+mult)*linofs) c CALL setviewport(1,linofs+1,maxx-1,maxy-mult*linofs-1) goto 50 endif C C...Label atoms, the floating point coordinates of atoms first have to be C converted to pixel coordinates, then the label is written using the C plotter character routine C $NODEBUG CALL setlinestyle(#FFFF) $DEBUG IF(Inumbr.lt.1)GOTO 50 DO 51 K=1,NATMS ISTART=maxx*(ixrev*CORPLT(K,IX)-XMIN)/RANGEX JSTART=IYSPAN-IYSPAN*(iyrev*CORPLT(K,IY)-YMIN)/RANGEY+2 IF(K.LT.10)THEN WRITE(LINE,'(I1)')K ISTART=ISTART-1 ENDIF IF(K.GE.10.AND.K.LT.100)THEN WRITE(LINE,'(I2)')K ISTART=ISTART-4 ENDIF IF(K.GE.100)THEN WRITE(LINE,'(I3)')K ISTART=ISTART-8 ENDIF IF(WT(K).GT.0.5D0)THEN dummy=setcolor(0) ELSE dummy=setcolor(15) ENDIF CALL PLOTCH(1,1,ISTART,JSTART,1) 51 CONTINUE dummy=setcolor(15) C 50 RETURN END C C_____________________________________________________________________________ C SUBROUTINE ROTATE(ikeepc) C INCLUDE 'FGRAPH.FD' MSF5+PS1 c USE MSFLIB PS4 c RECORD /rccoord/curpos RECORD /xycoord/ixy RECORD /wxycoord/wxy C logical*2 true PARAMETER (NATMAX=500,NIBOND=4*NATMAX,true=.true.) COMMON /limits/maxx,maxy,linofs,curpos,ixy,wxy, * mymode,myrows,mycols,blue,red COMMON /BLK1/NATMS,CORD,CORPLT,AXEND,ASPECT/BLK2/IBOND,NBONDS,WT, 1 radius/PARBLK/RANGEA,RANGEB,AMIN,BMIN COMMON /lines/toplin,botlin,emplin,COMENT REAL*8 CORD(NATMAX,3),CORPLT(NATMAX,3),WT(NATMAX),AXEND(6,3), * ASPECT,radius(natmax) REAL*8 ANGLE,RANGEA,RANGEB,SMALLA,SMALLB,AMIN,AMAX,BMIN,BMAX, 1 AMEAN,BMEAN INTEGER*2 maxx,maxy,mymode,myrows,mycols,linofs,dummy INTEGER*4 dummy4,blue,red INTEGER*2 IBOND(NIBOND,2),INKEY,N CHARACTER*79 toplin,botlin,emplin,COMENT CHARACTER KK,OUTSTR*21 C 1 WRITE(*,'(/'' Specify the projection plane IX,IY'',t37, * ''- use 1,2,3 for a,b,c''/t37, * ''- negative value reverses axis direction''// * t31,''..... '',$)') READ(*,*,ERR=1)IA,IB IF(iabs(IA).EQ.iabs(IB))GOTO 1 IF(iabs(IA).LT.1.OR.iabs(IA).GT.3)GOTO 1 IF(iabs(IB).LT.1.OR.iabs(IB).GT.3)GOTO 1 IC=6-(iabs(IA)+iabs(IB)) icontr=12 iarev=isign(1,ia) ibrev=isign(1,ib) icrev=1 ia=iabs(ia) ib=iabs(ib) C C---P-L-O-T---, first set up graphics conditions C C...Determine the graphics card and set the highest resolution graphics mode. C Use the second graphics page (page 1), this is convenient with Hercules C but may have to be changed to page 0 with some boards. Also determine C pixel colour for use in erasing points from previous pass C 333 call graphicsmode() C C...It is assumed that the ratio of physical screen dimensions Y/X is 0.75. C The pixel scaling ratio along the Y-coordinate for 1:1 X:Y dispersion C is maxy/(0.75*maxx) and is a unity for EGA and VGA, and 0.645 for C Hercules. C ASPECT=maxy/(0.75*maxx) C C...The setlinestyle parameter is INTEGER*2 and for continuous line has to C be 65535 - this is identified as an error by $DEBUG since I*2 range is C +-32767, hence have to $NODEBUG locally C $NODEBUG CALL setlinestyle(#FFFF) $DEBUG C C...Write top and bottom of display C if(myrows.eq.25.or.myrows.eq.30)then mult=1 else mult=2 endif CALL moveto(0,linofs,ixy) dummy=lineto(maxx,linofs) dummy=lineto(maxx,maxy-mult*linofs) dummy=lineto(0,maxy-mult*linofs) dummy=lineto(0,linofs) C WRITE(toplin,'(a)')coment WRITE(botlin,'(2A)')'Rotations: S-D E-X W-R ', 1 '|G|le |A|xes |N|umbers |P|arameters |H|elp' C 33 CALL settextposition(1,1,curpos) CALL outtext(toplin) CALL settextposition(myrows,1,curpos) CALL outtext(botlin) C C...Plot selected projection. The first step is scaling and since this code C is executed on each rotation the scaling varies continuously for maximum C screen filling C 3 CALL SCALIN(IA,iarev,RANGEA,SMALLA) CALL SCALIN(IB,ibrev,RANGEB,SMALLB) AMEAN=SMALLA+0.5D0*RANGEA BMEAN=SMALLB+0.5D0*RANGEB RANGEA=1.05D0*RANGEA RANGEB=1.05D0*RANGEB PIXELR=REAL(maxy-(1+mult)*linofs)/REAL(maxx) IF(RANGEB*ASPECT/RANGEA.LE.PIXELR)THEN RANGEB=RANGEA*PIXELR/ASPECT ELSE RANGEA=RANGEB*ASPECT/PIXELR ENDIF AMIN=AMEAN-RANGEA*0.5D0 AMAX=AMEAN+RANGEA*0.5D0 BMIN=BMEAN-0.5D0*RANGEB BMAX=BMEAN+0.5D0*RANGEB C C...Define graphics viewport and floating point window for automatic C coordinate to pixel conversion. Note that setviewport C puts the origin of the pixel axes at (0,linofs+1). C CALL setviewport(1,linofs+1,maxx-1,maxy-mult*linofs-1) dummy=setwindow(TRUE,amin,bmin,amax,bmax) C C...Plot molecule CALL VIEW(IA,IB,iarev,ibrev,ICONTR,1,ikeepc) C C---O-P-T-I-O-N-S--- C 2 IK=INKEY(N) IF(IK.EQ.13)GOTO 20 IF(IK.LT.65)GOTO 2 ANGLE=5. IF(IK.GT.90)THEN ANGLE=1. IK=IK-32 ENDIF KK=CHAR(IK) IF(KK.EQ.'E'.OR.KK.EQ.'X')GOTO 5 IF(KK.EQ.'S'.OR.KK.EQ.'D')GOTO 6 IF(KK.EQ.'W'.OR.KK.EQ.'R')GOTO 7 IF(KK.EQ.'A')then if(icontr-(icontr/10)*10.eq.0)then icontr=icontr+2 else icontr=icontr-2 endif CALL VIEW(IA,IB,iarev,ibrev,icontr,1,ikeepc) endif IF(KK.EQ.'N')then if(icontr/10.lt.1)then icontr=icontr+10 else icontr=icontr-10 endif CALL VIEW(IA,IB,iarev,ibrev,icontr,1,ikeepc) endif IF(KK.EQ.'G')then CALL VIEW(-IA,IB,iarev,ibrev,12,0,ikeepc) goto 333 endif IF(KK.EQ.'H')GOTO 500 IF(KK.EQ.'P')THEN if(icontr/10.lt.1)icontr=icontr+10 CALL VIEW(IA,IB,iarev,ibrev,icontr,1,ikeepc) CALL PARCHE GOTO 33 ENDIF GOTO 2 C C Exit C 20 dummy4 = setbkcolor( RED ) CALL settextposition(1,1,curpos) CALL outtext(emplin) WRITE(OUTSTR,'(A)')'ARE YOU SURE ?' CALL settextposition(1,1,curpos) CALL outtext(outstr(1:14)) CALL settextposition(1,20,curpos) WRITE(*,'(1X,A1,$)')CHAR(7) 916 IK=INKEY(N) KK=CHAR(IK) IF(KK.EQ.'Y'.OR.KK.EQ.'y')GOTO 915 IF(KK.NE.'N'.AND.KK.NE.'n')GOTO 916 dummy4 = setbkcolor( BLUE ) CALL settextposition(1,1,curpos) CALL outtext(toplin) GOTO 2 915 dummy=setvideomode($DEFAULTMODE) GOTO 4 C C Help screen C 500 dummy4 = setbkcolor( BLUE ) CALL clearscreen($GCLEARSCREEN) CALL settextposition(3,1,curpos) WRITE(*,232) 232 FORMAT( ' SUMMARY OF COMMANDS ACTIVE IN GRAPHICS MODE:'/ * 2x,44(1h-)// *' E/N S/D W/R - rotation of molecule about principal axes'/ *' (upper/lower case = fast/slow)'/ *' A - toggle plotting of principal axes'/ *' N - toggle numbering of atoms'/ *' G - write GLE copy of the diplayed orientation'/ *' H - display this help screen'// *' P - check structural parameters. When in this mode'/ *' B selects bond length, A bond angle and D -'/ *' dihedral angle. Once parameter type has been'/ *' selected numbers of atoms are to be typed in. ', * ' One digit'/ *' numbers are to be terminated with , two'/ *' digit numbers need no termination'// *' Exit from parameter check mode is with E, exit from graphics '/ *' mode is with followed by Y'/) WRITE(*,502) 502 FORMAT(1X,40X,' Press ENTER to return to diagram ',$) 501 IK=INKEY(N) IF(IK.EQ.0)GOTO 501 GOTO 333 C C Rotation about the horizontal axis C 5 IF(IK.EQ.69)ANGLE=-ANGLE CALL TRANSF(ANGLE,IB,IC,ibrev,icrev) GOTO 3 C C Rotation about the vertical axis C 6 IF(IK.EQ.83)ANGLE=-ANGLE CALL TRANSF(ANGLE,IA,IC,iarev,icrev) GOTO 3 C C Rotation about the axis perpendicular to the projection plane C 7 IF(IK.EQ.87)ANGLE=-ANGLE CALL TRANSF(ANGLE,IA,IB,iarev,ibrev) GOTO 3 C 4 RETURN END C C_____________________________________________________________________________ C SUBROUTINE PARCHE C C Check of structural parameters (bond lengths and bond angles) using C a scrollable window for listing the results C INCLUDE 'FGRAPH.FD' MSF5+PS1 c USE MSFLIB PS4 c RECORD /rccoord/curpos RECORD /xycoord/ixy RECORD /wxycoord/wxy C PARAMETER (NATMAX=500) COMMON /limits/maxx,maxy,linofs,curpos,ixy,wxy, * mymode,myrows,mycols,blue,red COMMON /BLK1/NATMS,CORD,CORPLT,AXEND,ASPECT COMMON /PARBLK/RANGEA,RANGEB,AMIN,BMIN COMMON /lines/toplin,botlin,emplin,COMENT CHARACTER KK,outlin*20,TEXT*20,RTEXT*21 REAL*8 BONDL,ANGLE,DANGLE,CORD(NATMAX,3),CORPLT(NATMAX,3), * AXEND(6,3) REAL*8 ANG,CON,COSANG,RA,RANGEA,RANGEB, 1 ASPECT,AMIN,BMIN INTEGER*2 maxx,maxy,mymode,myrows,mycols,linofs,inkey,dummy INTEGER*4 blue,red INTEGER*2 NPOS,NBUF CHARACTER*79 toplin,botlin,emplin,COMENT CHARACTER*36 RES(50),CURRES,BLANK C CON=0.174532925D-1 NBUF=0 rtext='+- scrollback |E|xit' WRITE(BLANK,'(36(1H ))') CALL settextposition(myrows,1,curpos) CALL outtext(emplin) C C C INPUT LOOP C C...selection of B (bond), A (angle), D (dihedral angle), C scrollback or E (exit) C 3 CALL settextposition(1,1,curpos) CALL outtext(emplin) CALL settextposition(1,1,curpos) CALL outtext('|B|ond,|A|ngle,|D|ihangle?') CALL settextposition(1,58,curpos) CALL outtext(rtext) 33 IK=INKEY(dummy) IF(IK.EQ.0)GOTO 33 IF(IK.GT.90)IK=IK-32 KK=CHAR(IK) IF(IK.EQ.'+'.or.KK.eq.'='.or.KK.eq.'-'.or.KK.eq.'_')GOTO 8 IF(KK.NE.'A'.AND.KK.NE.'B'.AND.KK.NE.'D'.AND.KK.NE.'E')GOTO 4 IF(KK.EQ.'E')GOTO 5 C CALL settextposition(1,1,curpos) CALL outtext(emplin) CALL settextposition(1,1,curpos) if(kk.eq.'B')CALL outtext( * '--> Type in numbers of atoms M,N for bond MN') if(kk.eq.'A')CALL outtext( * '--> Type in numbers of atoms or M,N,a for angle MN.a etc. ') if(kk.eq.'D')CALL outtext( * '--> Type in numbers of atoms K,L,M,N for dihedral angle KLMN') CALL settextposition(myrows,1,curpos) CALL outtext(BLANK) CALL settextposition(myrows,1,curpos) IF(KK.EQ.'B')THEN TEXT='----------Bond(' ELSE IF(KK.EQ.'A')THEN TEXT='---------Angle(' ELSE IF(KK.EQ.'D')THEN TEXT='Dihedral Angle(' ENDIF CALL outtext(TEXT) CURRES=TEXT NPOS=16 C C...read first two atom numbers, and calculate the relevant bond length CALL INNUM(N,NPOS,CURRES) CALL settextposition(myrows,NPOS,curpos) CALL outtext(',') CURRES=CURRES(1:(NPOS-1))//',' NPOS=NPOS+1 CALL INNUM(NN,NPOS,CURRES) IF(N.EQ.NN)GOTO 3 IF(N.LT.0.OR.NN.LT.0)GOTO 3 IF(KK.EQ.'B')THEN RA=BONDL(CORD(N,1),CORD(N,2),CORD(N,3), * CORD(NN,1),CORD(NN,2),CORD(NN,3)) WRITE(outlin,3104)RA 3104 FORMAT(') =',F8.5,$) CALL settextposition(myrows,NPOS,curpos) CALL outtext(outlin) GOTO 6 ENDIF C C...third atom number for angle and dihedral angle C (if angle selected then a,b or c are also allowed on input for C determining angles between bonds and principal axes) C CALL settextposition(myrows,NPOS,curpos) CALL outtext(',') CURRES=CURRES(1:(NPOS-1))//',' NPOS=NPOS+1 CALL INNUM(NNN,NPOS,CURRES) IF(NNN.EQ.N.OR.NNN.EQ.NN)GOTO 3 IF(NNN.LT.0.AND.KK.NE.'A')GOTO 3 IF(KK.EQ.'A')THEN C C...angle between bond and axis IF(NNN.LT.0)THEN NNN=-NNN RA=BONDL(CORD(N,1),CORD(N,2),CORD(N,3), * CORD(NN,1),CORD(NN,2),CORD(NN,3)) COSANG=(CORD(NN,NNN)-CORD(N,NNN))/RA ANG=DACOS(COSANG)/CON C C...angle between two bonds ELSE ANG=ANGLE(CORD(N,1),CORD(N,2),CORD(N,3), * CORD(NN,1),CORD(NN,2),CORD(NN,3), * CORD(NNN,1),CORD(NNN,2),CORD(NNN,3)) ENDIF C WRITE(outlin,3103)ANG 3103 FORMAT(') =',F8.3,$) CALL settextposition(myrows,NPOS,curpos) CALL outtext(outlin) GOTO 6 ENDIF C C...fourth atom number for dihedral angle CALL settextposition(myrows,NPOS,curpos) CALL outtext(',') CURRES=CURRES(1:(NPOS-1))//',' NPOS=NPOS+1 CALL INNUM(NNNN,NPOS,CURRES) IF(NNNN.EQ.N.OR.NNNN.EQ.NN.OR.NNNN.EQ.NNN)GOTO 3 IF(NNNN.LT.0)GOTO 3 ANG=DANGLE(CORD(N,1),CORD(N,2),CORD(N,3), * CORD(NN,1),CORD(NN,2),CORD(NN,3), * CORD(NNN,1),CORD(NNN,2),CORD(NNN,3), * CORD(NNNN,1),CORD(NNNN,2),CORD(NNNN,3)) WRITE(outlin,3103)ANG CALL settextposition(myrows,NPOS,curpos) CALL outtext(outlin) GOTO 6 C C...Beep for input error and return to option selection C 4 CALL settextposition(1,1,curpos) WRITE(*,'(1X,A1,$)')CHAR(7) GOTO 3 C C...Update the scrollback buffer C 6 CURRES=CURRES(1:(NPOS-1))//outlin DO 7 N=19,1,-1 RES(N+1)=RES(N) 7 CONTINUE RES(1)=CURRES NBUF=NBUF+1 GOTO 3 C C...Scroll through the results buffer C 8 IF(KK.EQ.'+'.or.KK.eq.'=')THEN NBUF=NBUF+1 IF(NBUF.GT.50)NBUF=1 ELSE NBUF=NBUF-1 IF(NBUF.LT.1)NBUF=50 ENDIF CALL settextposition(myrows,43,curpos) CALL outtext(BLANK) CALL settextposition(myrows,43,curpos) CALL outtext(RES(NBUF)) GOTO 3 C 5 RETURN END C C_____________________________________________________________________________ C SUBROUTINE PLOTCH(XSCALE,YSCALE,ISTART,JSTART,IDIR) C C This is a simple plotter type character generator which plots the C specified character string LINE starting from absolute pixel coordinates C (ISTART,JSTART) The font is separately scalable in horizontal and C vertical directions and choice of upright type or italics is possible. C C LINE - text string to be written C XSCALE - integer scaling factor for the horizontal character dimension C YSCALE - integer scaling factor for the vertical character dimension C ISTART, JSTART - starting coordinates for placing the text in units of C the output device C IDIR - direction of writing the character string C = +-1 left to right, = +-2 from bottom to top C (negative values specify italics) C C INCLUDE 'FGRAPH.FD' MSF5+PS1 c USE MSFLIB PS4 c RECORD /rccoord/curpos RECORD /xycoord/ixy RECORD /wxycoord/wxy C COMMON /limits/maxx,maxylinofs,curpos,ixy,wxy, * mymode,myrows,mycols,blue,red COMMON /TXTBUF/LINE EQUIVALENCE (LINE,TEXT(1)) EQUIVALENCE(CHARPL(1,1),ACHAR(1,1)),(CHARPL(1,11),BCHAR(1,1)), 1 (CHARPL(1,21),CCHAR(1,1)),(CHARPL(1,31),DCHAR(1,1)), 1 (CHARPL(1,41),ECHAR(1,1)),(CHARPL(1,51),FCHAR(1,1)), 1 (CHARPL(1,61),GCHAR(1,1)),(CHARPL(1,71),HCHAR(1,1)), 1 (CHARPL(1,81),JCHAR(1,1)),(CHARPL(1,91),KCHAR(1,1)) CHARACTER*80 LINE CHARACTER TEXT(80),NULsym,SPACE INTEGER*2 maxx,maxy,mymode,myrows,mycols,linofs,dummy INTEGER*4 blue,red INTEGER*4 XSCALE,YSCALE,ISTART,JSTART,IDIR INTEGER*2 NPLOT(95),IX,IY,IXX,IYY INTEGER*2 ACHAR(22,10),BCHAR(22,10),CCHAR(22,10),DCHAR(22,10), 1 ECHAR(22,10),FCHAR(22,10),GCHAR(22,10),HCHAR(22,10), 1 JCHAR(22,10),KCHAR(22,5),CHARPL(22,95) C SLANT=0. IDIREC=IABS(IDIR) IXSTRT=ISTART IYSTRT=JSTART WRITE(NULsym,'(A1)')CHAR(0) SPACE=' ' DO 100 NCHARS=80,1,-1 100 IF(TEXT(NCHARS).NE.NULsym.AND.TEXT(NCHARS).NE.SPACE)GOTO 101 C 101 GO TO (31,32),IDIREC 31 SCALEX=XSCALE SCALEY=YSCALE ISPACE=3*XSCALE IF(IDIR.LT.0)SLANT=0.5*SCALEY GO TO 30 32 SCALEX=YSCALE*1.5 SCALEY=XSCALE IF(SCALEY.GT.1.)SCALEY=SCALEY/1.5 ISPACE=2*XSCALE IF(IDIR.LT.0)SLANT=0.334*SCALEY 30 CONTINUE C C CHARACTER DEFINITIONS - NPLOT contains the number of points in character C definition, CHARPL contains the definitions C C ASCII codes from 32 to 126 are supported with their usual meanings C with the following exceptions: C ^ ( 94) - shift up for superscripts C ~ (126) - shift down for subscripts C ` ( 96) - one half (1 over 2) character C @ ( 64) - arrow to the left ie. <- C C _!"#$%&'() *+,-./0123 456789:;<= >?@ABCDEFG HIJKLMNOPQ C C RSTUVWXYZ[ \]^_`abcde fghijklmno pqrstuvwxy z{|}~ C DATA NPLOT/6,4,6,8,10,1,11,6,4,4, * 10,4,6,2,5,2,9,5,7,9, 4 5,8,11,4,11,11,10,11,3,4, > 3,10,10,5,11,8,7,6,5,10, H 6,6,5,7,3,5,4,9,7,11, R 9,10,4,6,3,5,4,5,4,4, \ 2,4,1,2,6,9,10,8,11,9, f 6,10,7,7,6,7,5,10,7,9, p 10,10,6,8,6,8,3,11,4,6, z 4,9,2,9,1/ C DATA ACHAR/0,0,1,6,3,6,4,5,4,4,3,3,0,0,0,0,0,0,0,0,0,0, space ! 0,6,0,2,100,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, " 1,6,1,5,0,4,103,6,3,5,2,4,0,0,0,0,0,0,0,0,0,0, # 1,6,1,0,103,6,3,0,100,4,4,4,100,2,4,2,0,0,0,0,0,0, $ 2,6,2,0,100,1,3,1,4,2,3,3,1,3,0,4,1,5,4,5,0,0, % 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, & 4,0,0,4,0,5,1,6,2,5,2,4,0,2,0,1,1,0,2,0,4,2, ' 1,5,0,5,0,6,1,6,1,4,0,3,0,0,0,0,0,0,0,0,0,0, ( 2,7,0,5,0,1,2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0, ) 0,7,2,5,2,1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ DATA BCHAR/1,5,2,4,2,2,1,1,103,5,2,4,102,2,3,1,100,3,4,3,0,0, * + 0,3,4,3,102,5,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0, , 1,0,0,0,0,1,1,1,1,-1,0,-2,0,0,0,0,0,0,0,0,0,0, - 0,3,4,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, . 0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, / 0,0,4,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, O 1,0,3,0,4,1,4,5,3,6,1,6,0,5,0,1,1,0,0,0,0,0, 1 0,0,2,0,101,0,1,6,0,5,0,0,0,0,0,0,0,0,0,0,0,0, 2 4,0,0,0,4,4,4,5,3,6,1,6,0,5,0,0,0,0,0,0,0,0, 3 0,6,4,6,2,4,3,4,4,3,4,1,3,0,1,0,0,1,0,0,0,0/ DATA CCHAR/3,0,3,6,0,3,0,2,4,2,0,0,0,0,0,0,0,0,0,0,0,0, 4 5 4,6,0,6,0,3,3,3,4,2,4,1,3,0,0,0,0,0,0,0,0,0, 6 3,6,2,6,0,4,0,1,1,0,3,0,4,1,4,2,3,3,1,3,0,2, 7 0,6,4,6,1,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 8 1,3,3,3,4,2,4,1,3,0,1,0,0,1,0,2,1,3,0,4,0,5, 9 1,0,2,0,4,2,4,5,3,6,1,6,0,5,0,4,1,3,3,3,4,4, : 0,4,0,5,1,5,1,4,0,4,100,0,1,0,1,1,0,1,0,0,0,0, ; 0,4,0,5,1,5,1,4,0,4,101,0,0,0,0,1,1,1,1,-1,0,-2, < 3,6,0,3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, = 0,4,4,4,100,2,4,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ DATA DCHAR/0,6,3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, > ? 0,5,1,6,3,6,4,5,4,4,3,3,2,3,2,2,102,0,2,0,0,0, @ 0,7,1,8,1,5,100,3,3,3,100,1,1,1,2,0,0,-2,2,-2,0,0, A 0,0,2,6,4,0,101,2,3,2,0,0,0,0,0,0,0,0,0,0,0,0, B 0,3,3,3,4,4,4,5,3,6,0,6,0,0,3,0,4,1,4,2,3,3, C 4,5,3,6,1,6,0,5,0,1,1,0,3,0,4,1,0,0,0,0,0,0, D 0,0,3,0,4,1,4,5,3,6,0,6,0,0,0,0,0,0,0,0,0,0, E 4,0,0,0,0,6,4,6,100,3,2,3,0,0,0,0,0,0,0,0,0,0, F 0,0,0,6,4,6,100,3,2,3,0,0,0,0,0,0,0,0,0,0,0,0, G 4,5,3,6,1,6,0,5,0,1,1,0,3,0,4,1,4,3,2,3,0,0/ DATA ECHAR/0,0,0,6,100,3,4,3,104,6,4,0,0,0,0,0,0,0,0,0,0,0, H I 0,0,2,0,101,0,1,6,100,6,2,6,0,0,0,0,0,0,0,0,0,0, J 0,0,2,0,3,1,3,6,0,6,0,0,0,0,0,0,0,0,0,0,0,0, K 0,0,0,6,100,2,2,3,4,6,102,3,4,0,0,0,0,0,0,0,0,0, L 0,6,0,0,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, M 0,0,0,6,3,3,6,6,6,0,0,0,0,0,0,0,0,0,0,0,0,0, N 0,0,0,6,4,0,4,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0, O 1,0,3,0,4,1,4,5,3,6,1,6,0,5,0,1,1,0,0,0,0,0, P 0,0,0,6,3,6,4,5,4,3,3,2,0,2,0,0,0,0,0,0,0,0, Q 1,0,3,0,4,1,4,5,3,6,1,6,0,5,0,1,1,0,103,1,4,-1/ DATA FCHAR/0,0,0,6,3,6,4,5,4,3,3,2,0,2,102,2,4,0,0,0,0,0, R S 4,6,1,6,0,5,0,4,1,3,3,3,4,2,4,1,3,0,0,0,0,0, T 2,0,2,6,100,6,4,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0, U 0,6,0,1,1,0,3,0,4,1,4,6,0,0,0,0,0,0,0,0,0,0, V 0,6,2,0,4,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, W 0,6,0,0,3,3,6,0,6,6,0,0,0,0,0,0,0,0,0,0,0,0, X 0,6,4,0,100,0,4,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0, Y 2,0,2,3,0,6,104,6,2,3,0,0,0,0,0,0,0,0,0,0,0,0, Z 4,0,0,0,4,6,0,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0, [ 2,7,0,7,0,-1,2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ DATA GCHAR/0,6,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, \ ] 0,-1,2,-1,2,7,0,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0, ^ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, _ 0,-2,4,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, ` 0,5,1,5,1,6,0,6,0,4,1,3,0,0,0,0,0,0,0,0,0,0, a 4,1,3,0,1,0,0,1,0,3,1,4,4,4,4,1,5,0,0,0,0,0, b 0,6,0,0,100,1,1,0,3,0,4,1,4,3,3,4,1,4,0,3,0,0, c 4,1,3,0,1,0,0,1,0,3,1,4,3,4,4,3,0,0,0,0,0,0, d 4,6,4,1,5,0,104,1,3,0,1,0,0,1,0,3,1,4,3,4,4,3, e 3,0,1,0,0,1,0,3,1,4,3,4,4,3,4,2,0,2,0,0,0,0/ DATA HCHAR/1,0,1,5,2,6,3,5,100,3,2,3,0,0,0,0,0,0,0,0,0,0, f g 0,-2,3,-2,4,-1,4,4,1,4,0,3,0,1,1,0,3,0,4,1,0,0, h 0,0,0,6,100,3,1,4,3,4,4,3,4,0,0,0,0,0,0,0,0,0, i 0,4,1,4,1,0,100,0,2,0,101,6,1,6,0,0,0,0,0,0,0,0, j 0,-2,1,-2,2,-1,2,4,102,6,2,6,0,0,0,0,0,0,0,0,0,0, k 0,0,0,6,100,1,2,2,4,4,102,2,4,0,0,0,0,0,0,0,0,0, l 0,6,1,6,1,0,100,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0, m 0,0,0,4,100,3,2,4,3,4,3,0,103,3,5,4,6,4,6,0,0,0, n 0,0,0,4,100,3,2,4,3,4,4,3,4,0,0,0,0,0,0,0,0,0, o 1,0,3,0,4,1,4,3,3,4,1,4,0,3,0,1,1,0,0,0,0,0/ DATA JCHAR/0,-2,0,4,100,1,1,0,3,0,4,1,4,3,3,4,1,4,0,3,0,0, p q 4,-2,4,4,104,3,3,4,1,4,0,3,0,1,1,0,3,0,4,1,0,0, r 0,0,0,4,100,2,2,4,3,4,4,3,0,0,0,0,0,0,0,0,0,0, s 4,4,1,4,0,3,1,2,3,2,4,1,3,0,0,0,0,0,0,0,0,0, t 1,6,1,1,2,0,3,0,100,4,2,4,0,0,0,0,0,0,0,0,0,0, u 0,4,0,1,1,0,3,0,4,1,4,4,104,1,5,0,0,0,0,0,0,0, v 0,4,2,0,4,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, w 0,4,0,1,1,0,2,0,3,1,3,4,103,1,4,0,5,0,6,1,6,4, x 0,0,4,4,100,4,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, y 0,4,2,0,104,4,2,0,1,-2,0,-2,0,0,0,0,0,0,0,0,0,0/ DATA KCHAR/0,4,4,4,0,0,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, z { 3,7,2,7,1,6,1,4,0,3,1,2,1,0,2,-1,3,-1,0,0,0,0, | 0,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, } 0,-1,1,-1,2,0,2,2,3,3,2,4,2,6,1,7,0,7,0,0,0,0, ~ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ C C C character identification and writing loop C ILIFT=0 DO 3 N=1,NCHARS IK=ICHAR(TEXT(N)) IF(IK.GT.32.AND.IK.LE.126)GOTO 6 IF(IK.NE.32)GO TO 3 C C...... blank character GO TO (33,34),IDIREC 33 IXSTRT=IXSTRT+5*SCALEX GO TO 3 34 IYSTRT=IYSTRT-5*SCALEY GO TO 3 C 6 NN=IK-31 IF(IK.NE.126.AND.IK.NE.94)GOTO 5 C C...... subscripts and superscripts ISHIFT=4 IF(IK.EQ.126)ISHIFT=-4 GOTO (115,116),IDIREC 115 IYSTRT=IYSTRT-ISHIFT*SCALEY IXSTRT=IXSTRT+ISPACE*SCALEX*ISHIFT/16 GOTO 3 116 IXSTRT=IXSTRT-ISHIFT*SCALEX IYSTRT=IYSTRT-ISPACE*SCALEY*ISHIFT/16 GOTO 3 C 5 NPTS=NPLOT(NN) C C..... move to beginning of character C GO TO (35,36),IDIREC 35 IX=IXSTRT+CHARPL(1,NN)*SCALEX+CHARPL(2,NN)*SLANT IY=IYSTRT-CHARPL(2,NN)*SCALEY GO TO 37 36 IX=IXSTRT-CHARPL(2,NN)*SCALEX IY=IYSTRT-CHARPL(1,NN)*SCALEY-CHARPL(2,NN)*SLANT 37 call moveto(ix,iy,ixy) C C...... go through points defining the character C 114 DO 7 MN=3,2*NPTS,2 I=CHARPL(MN,NN) J=CHARPL(MN+1,NN) IF(I.LT.90)GOTO 110 ILIFT=1 I=I-100 110 GO TO (38,39),IDIREC 38 IXX=IXSTRT+I*SCALEX+J*SLANT IYY=IYSTRT-J*SCALEY GO TO 40 39 IXX=IXSTRT-J*SCALEX IYY=IYSTRT-I*SCALEY-J*SLANT C C...... move pen to next point of character definition C 40 IF(ILIFT.NE.1)GOTO 111 call moveto(IXX,IYY,ixy) ILIFT=0 GOTO 7 111 dummy=lineto(IXX,IYY) C 7 CONTINUE C C Backtrack for remainder of character 8 (16 defining points) C IF(IK.NE.56)GOTO 112 IK=32 NN=1 NPTS=6 GOTO 114 C C Find starting position for next character C 112 LARGE=0 DO 9 MMNN=1,2*NPTS,2 J=CHARPL(MMNN,NN) IF(J.GT.90)J=J-100 9 IF(J.GT.LARGE)LARGE=J GO TO(41,42),IDIREC 41 IXSTRT=IXSTRT+ISPACE+LARGE*SCALEX GO TO 3 42 IYSTRT=IYSTRT-ISPACE-LARGE*SCALEY C 3 CONTINUE C C RETURN END C C------------------------------------------------------------------------ C SUBROUTINE INNUM(N,NPOS,CURRES) C C Input and check of validity of a 1 to 2 digit integer number; C N - takes on the value of the input atom number (or -1,-2,-3 for C a,b,c axis in the case of angle between a bond and a principal C axis) C NPOS - starting column for writing of next number to bottom line of C screen, updated by INNUM after the number is written C CURRES - buffer containing the current result line C INCLUDE 'FGRAPH.FD' MSF5+PS1 c USE MSFLIB PS4 c RECORD /rccoord/curpos RECORD /xycoord/ixy RECORD /wxycoord/wxy C PARAMETER (NATMAX=500) COMMON /limits/maxx,maxy,linofs,curpos,ixy,wxy, * mymode,myrows,mycols,blue,red COMMON /BLK1/NATMS,COR,CORPLT,AXEND,ASPECT CHARACTER TEXT*2,KK,CURRES*36 INTEGER*2 maxx,maxy,linofs,mymode,myrows,mycols,dummy,inkey, * NDIG,IK,NPOS INTEGER*4 blue,red REAL*8 COR(NATMAX,3),CORPLT(NATMAX,3),AXEND(6,3),ASPECT NDIG=0 C C...Accept input C CALL settextposition(myrows,NPOS,curpos) 5 IK=inkey(dummy) IF(IK.EQ.0)GOTO 5 KK=CHAR(IK) C C...identify axis IF(KK.EQ.'A'.OR.KK.EQ.'a')THEN CALL outtext('a') CURRES=CURRES(1:(NPOS-1))//'a' NPOS=NPOS+1 N=-1 GOTO 4 ENDIF IF(KK.EQ.'B'.OR.KK.EQ.'b')THEN CALL outtext('b') CURRES=CURRES(1:(NPOS-1))//'b' NPOS=NPOS+1 N=-2 GOTO 4 ENDIF IF(KK.EQ.'C'.OR.KK.EQ.'c')THEN CALL outtext('c') CURRES=CURRES(1:(NPOS-1))//'c' NPOS=NPOS+1 N=-3 GOTO 4 ENDIF C C...identify number (one or two digit) IF(IK.EQ.13.AND.NDIG.EQ.1)THEN N=N1 KK=CHAR(N+48) IF(N.GT.NATMS)THEN WRITE(*,'(1X,A1,$)')CHAR(7) NDIG=0 GOTO 5 ENDIF CALL outtext(KK) CURRES=CURRES(1:(NPOS-1))//KK NPOS=NPOS+1 GOTO 4 ENDIF IF(KK.LT.'0'.OR.KK.GT.'9')GOTO 5 IF(NDIG.EQ.0)THEN N1=IK-48 NDIG=1 GOTO 5 ENDIF IF(NDIG.EQ.1)THEN NDIG=2 N=N1*10+(IK-48) IF(N.GT.NATMS)THEN WRITE(*,'(1X,A1,$)')CHAR(7) NDIG=0 GOTO 5 ENDIF WRITE(TEXT,'(I2)')N CALL outtext(TEXT) CURRES=CURRES(1:(NPOS-1))//TEXT NPOS=NPOS+2 GOTO 4 ENDIF C 4 RETURN END C C_____________________________________________________________________________ C SUBROUTINE NFIRST(LINE,IWHICH,REALNM,INTNUM,IERROR) C C This routine reads the first (integer or real) number from C the buffer and eliminates it from the buffer C C IWHICH=1 first integer number to be extracted C IWHICH=2 first real number to be extracted C PARAMETER (LINLEN=79) REAL*8 REALNM CHARACTER LINE*79,INFORM*7 C IERROR=0 INUM=0 DO 1 I=1,LINLEN IF((LINE(I:I).EQ.' '.OR.LINE(I:I).EQ.CHAR(0)).AND.INUM.EQ.1) * GOTO 2 IF(LINE(I:I).NE.' ')INUM=1 1 CONTINUE C 3 WRITE(*,'(1X/A/1X,A)')' **** INPUT PROBLEM ON PROCESSING:',LINE IERROR=1 RETURN C 2 IF(IWHICH.EQ.1)THEN WRITE(INFORM,'(''(I'',I2,'' )'')')I-1 READ(LINE,INFORM)INTNUM ENDIF IF(IWHICH.EQ.2)THEN WRITE(INFORM,'(''(F'',I2,''.0)'')')I-1 READ(LINE,INFORM)REALNM ENDIF LINE=LINE(I:) C RETURN END C C_____________________________________________________________________________ C SUBROUTINE CONVER C C Conversion from CART definitions to initial Cartesian coordinates C IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NATMAX=500) COMMON /BLK0/INITIA,NATOMS,ATOMS,VBA,TRANS,PRC,RJA,VCA 1 /BLK1/IFINAL,CORD,CORPLT,AXEND,ASPECT REAL*8 CORD(NATMAX,3),ATOMS(NATMAX,3),VBA(3),TRANS(3,3) 1 ,PRC(3),RJA(3),VCA(3),AXEND(6,3),CORPLT(NATMAX,3),ASPECT INTEGER*2 NATOMS(NATMAX,4) CON=0.174532925D-1 C DO 5556 I=INITIA,IFINAL NO=NATOMS(I,1) NA=NATOMS(I,2) NB=NATOMS(I,3) NC=NATOMS(I,4) IF(ATOMS(I,2))133,132,133 132 CS=-0.3333333333D0 SS=0.942809041D0 GO TO 5558 133 CS=DCOS(CON*ATOMS(I,2)) SS=DSIN(CON*ATOMS(I,2)) 5558 DSQ=0.0D0 DO 138 M=1,3 VBA(M)=CORD(NB,M)-CORD(NA,M) VCA(M)=CORD(NC,M)-CORD(NA,M) DSQ=DSQ+VBA(M)**2 138 CONTINUE RAB=DSQRT(DSQ) SCALE=0.0D0 DO 142 M=1,3 TRANS(M,1)=VBA(M)/RAB SCALE=SCALE+TRANS(M,1)*VCA(M) 142 CONTINUE DSQ=0.0D0 DO 146 M=1,3 RJA(M)=VCA(M)-SCALE*TRANS(M,1) DSQ=DSQ+RJA(M)**2 146 CONTINUE RAJ=DSQRT(DSQ) DO 148 M=1,3 TRANS(M,2)=RJA(M)/RAJ 148 CONTINUE TRANS(1,3)=TRANS(2,1)*TRANS(3,2)-TRANS(3,1)*TRANS(2,2) TRANS(2,3)=TRANS(3,1)*TRANS(1,2)-TRANS(1,1)*TRANS(3,2) TRANS(3,3)=TRANS(1,1)*TRANS(2,2)-TRANS(2,1)*TRANS(1,2) PRC(1)=ATOMS(I,1)*CS PRC(2)=ATOMS(I,1)*SS*DCOS(CON*ATOMS(I,3)) PRC(3)=ATOMS(I,1)*SS*DSIN(CON*ATOMS(I,3)) DO 5556 M=1,3 CORD(NO,M)=CORD(NA,M) DO 5556 K=1,3 CORD(NO,M)=CORD(NO,M)+TRANS(M,K)*PRC(K) 5556 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE BONLOC C C BONLOC locates bonds in the molecule ie. atoms that are closer together C than the sum of their covalent radii + 7% C C - numbers of atoms involved in an identified bond are placed in IBOND, C NBONDS is the number of bonds identified C - BONLOC will only recognise atoms of H,B,C,N,O,F,Si,P,S,Cl,Ge,As, C Se,Br,Sn,Sb,Te, and I in their stable isotopic species, ie. those C elements for which radii have been tabulated by Pauling. Some C minor isotopic species from Se onwards have been omitted where their C mass number clashes with more abundant isotope of another element. IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NRADII=45,NATMAX=500,NIBOND=4*NATMAX) COMMON /BLK1/NATMS,CORD,CORPLT,AXEND,ASPECT * /BLK2/IBOND,NBONDS,WT,radius REAL*8 CORD(NATMAX,3),WT(NATMAX),CORPLT(NATMAX,3),AXEND(6,3), * ASPECT,radius(natmax) INTEGER*2 IBOND(NIBOND,2),MASS(NRADII) REAL*4 CORAD(NRADII) C DATA MASS/ 1, 2, 10, 11, 12, 13, 14, 15, 16, 17, * 18, 19, 28, 29, 30, 31, 32, 33, 34, 35, 37, * 70, 72, 73, 74, 76, 75, 77, 78, 80, 82, 79, 81, * 116, 117, 118, 119, 120, 121, 123, 125, 126, 128, 130, 127/ DATA CORAD/0.37,0.37,0.81,0.81,0.77,0.77,0.74,0.74,0.71,0.71, * 0.71,0.72,1.17,1.17,1.17,1.10,1.04,1.04,1.04,0.99,0.99, * 1.22,1.22,1.22,1.22,1.22,1.21,1.17,1.17,1.17,1.17,1.14,1.14, * 1.40,1.40,1.40,1.40,1.40,1.41,1.41,1.37,1.37,1.37,1.37,1.33/ C 10 FORMAT(1X,'MASS OF ATOM.',I2,' = ',F8.4,' IS NOT INCLUDED', 1 ' IN BONLOC RECOGNITION TABLE') C NBONDS=0 DO 9 NSTART=2,NATMS NNN=NSTART-1 MMASS=ANINT(WT(NNN)) C C...treatment of dummy atoms IF(MMASS.EQ.0)THEN NBONDS=NBONDS+1 IBOND(NBONDS,1)=NNN IBOND(NBONDS,2)=NNN radius(nnn)=0.1d0 GOTO 9 ENDIF C C...treatment of normal atom pairs, first identify the current atom, then C go through remaining atoms and determine whether any are closer than sum C of covalent radii DO 3 NN=1,NRADII IF(MMASS.EQ.MASS(NN))GO TO 4 3 CONTINUE WRITE(*, 10)NNN,WT(NNN) NBONDS=NBONDS+1 IBOND(NBONDS,1)=NNN IBOND(NBONDS,2)=NNN GOTO 9 C 4 RAD=CORAD(NN) radius(nnn)=rad DO 1 N=NSTART,NATMS MMASS=IDINT(WT(N)/0.5)*0.5+0.5 IF(MMASS.EQ.0)GOTO 1 DO 2 NN=1,NRADII IF(MMASS.EQ.MASS(NN))GO TO 5 2 CONTINUE WRITE(*, 10)N,WT(N) GOTO 1 5 radius(n)=corad(nn) DIST=RAD+CORAD(NN) SEP=SQRT((CORD(NNN,1)-CORD(N,1))**2+(CORD(NNN,2)-CORD(N,2)) 1 **2 +(CORD(NNN,3)-CORD(N,3))**2) IF(SEP.LT.DIST*1.07)GO TO 6 GO TO 1 6 NBONDS=NBONDS+1 IBOND(NBONDS,1)=NNN IBOND(NBONDS,2)=N 1 CONTINUE 9 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE SCALIN(NAXIS,naxrev,RANGEX,SMALLX) C C Determination of minimum cooordinate value SMALLX and the range of c coordinate values RANGEX along a given principal axis NAXIS c NAXREV defines whether the axis is reversed (-1) or not (1) relative c to signs from diagonalization C C IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NATMAX=500) COMMON /BLK1/NATMS,CORD,CORPLT,AXEND,ASPECT REAL*8 CORD(NATMAX,3),CORPLT(NATMAX,3),AXEND(6,3),ASPECT C BIGX=-10.0 SMALLX=10.0 C DO 1 N=1,NATMS coord=naxrev*corplt(n,iabs(naxis)) IF(coord.GT.BIGX)BIGX=coord IF(coord.LT.SMALLX)SMALLX=coord 1 CONTINUE C RANGEX=1.1d0*(BIGX-SMALLX) smallx=1.05d0*smallx C RETURN END C_____________________________________________________________________________ C SUBROUTINE TRANSF(ANGLE,IX,IY,ixrev,iyrev) C C Transformation of coordinates IX,IY by rotation through an angle C of ANGLE degrees C IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NATMAX=500) COMMON /BLK1/NATMS,CORD,CORPLT,AXEND,ASPECT REAL*8 CORD(NATMAX,3),CORPLT(NATMAX,3),AXEND(6,3),ASPECT C CON=0.174532925D-1 ANGLE=ANGLE*CON S=DSIN(ANGLE) C=DCOS(ANGLE) C DO 1 N=1,NATMS A=C*ixrev*CORPLT(N,IX)-S*iyrev*CORPLT(N,IY) B=S*ixrev*CORPLT(N,IX)+C*iyrev*CORPLT(N,IY) CORPLT(N,IX)=ixrev*A CORPLT(N,IY)=iyrev*B 1 CONTINUE DO 2 N=1,6 A=C*ixrev*AXEND(N,IX)-S*iyrev*AXEND(N,IY) B=S*ixrev*AXEND(N,IX)+C*iyrev*AXEND(N,IY) AXEND(N,IX)=ixrev*A AXEND(N,IY)=iyrev*B 2 CONTINUE C RETURN END C C_____________________________________________________________________________ C FUNCTION BONDL(A1,A2,A3,B1,B2,B3) C C Bond (A.B) C REAL*8 A1,A2,A3,B1,B2,B3,BONDL C BONDL=DSQRT((A1-B1)**2+(A2-B2)**2+(A3-B3)**2) C RETURN END C_____________________________________________________________________________ C FUNCTION ANGLE(A1,A2,A3,B1,B2,B3,C1,C2,C3) C C Angle (A.B.C) using cosine rule, returned in degrees C REAL*8 A1,A2,A3,B1,B2,B3,C1,C2,C3,CONV,ANGLE,BONDL REAL*8 BL1,BL2,BL3 CONV=1.745329252D-2 C BL1=BONDL(A1,A2,A3,B1,B2,B3) BL2=BONDL(B1,B2,B3,C1,C2,C3) BL3=BONDL(A1,A2,A3,C1,C2,C3) ANGLE=(BL1**2+BL2**2-BL3**2)/(2.D0*BL1*BL2) IF(DABS(ANGLE).GT.1.D0)ANGLE=DSIGN(1.D0,ANGLE) ANGLE=DACOS(ANGLE)/CONV C RETURN END C_____________________________________________________________________________ C FUNCTION DANGLE(A1,A2,A3,B1,B2,B3,C1,C2,C3,D1,D2,D3) C C Dihedral angle (A.B.C.D), returned in degrees C REAL*8 A1,A2,A3,B1,B2,B3,C1,C2,C3,D1,D2,D3, * CONV,DANGLE,ANGLE,BONDL REAL*8 ANGL1,ANGL2,P1,P2,P3,Q1,Q2,Q3,R12S,R23, * R34S,BL1,BL2,DOTP CONV=1.745329252D-2 C ANGL1=ANGLE(A1,A2,A3,B1,B2,B3,C1,C2,C3)-90.D0 ANGL2=ANGLE(B1,B2,B3,C1,C2,C3,D1,D2,D3)-90.D0 R12S=BONDL(A1,A2,A3,B1,B2,B3)*DSIN(ANGL1*CONV) R23=BONDL(B1,B2,B3,C1,C2,C3) R34S=BONDL(C1,C2,C3,D1,D2,D3)*DSIN(ANGL2*CONV) C P1=A1-B1+R12S*(C1-B1)/R23 P2=A2-B2+R12S*(C2-B2)/R23 P3=A3-B3+R12S*(C3-B3)/R23 Q1=D1-C1-R34S*(C1-B1)/R23 Q2=D2-C2-R34S*(C2-B2)/R23 Q3=D3-C3-R34S*(C3-B3)/R23 C BL1=DSQRT(P1**2+P2**2+P3**2) BL2=DSQRT(Q1**2+Q2**2+Q3**2) DOTP=P1*Q1+P2*Q2+P3*Q3 C DANGLE=DOTP/(BL1*BL2) IF(DABS(DANGLE).GT.1.D0)DANGLE=DSIGN(1.D0,DANGLE) DANGLE=DACOS(DANGLE)/CONV C RETURN END C_____________________________________________________________________________ C SUBROUTINE EIGEN(A,R,N,MV) C C MATRIX DIAGONALIZATION ROUTINE FOR SYMMETRIC MATRICES FROM THE C 'SSP' LIBRARY C C - Jacobi type algorithm C - matrix to be diagonalized has to be in symmetric storage mode C (A and R are presently dimensioned for 10x10 matrices) C REAL*8 A(55),R(100) REAL*8 ANORM,ANRMX,THR,X,Y,SINX,SINX2,COSX, 1 COSX2,SINCS,RANGE,DSQRT,DABS C C C GENERATE IDENTITY MATRIX C 5 RANGE=1.0D-12 IF(MV-1) 10,25,10 10 IQ=-N DO 20 J=1,N IQ=IQ+N DO 20 I=1,N IJ=IQ+I R(IJ)=0.0 IF(I-J) 20,15,20 15 R(IJ)=1.0 20 CONTINUE C C COMPUTE INITIAL AND FINAL NORMS (ANORM AND ANORMX) C 25 ANORM=0.0 DO 35 I=1,N DO 35 J=I,N IF(I-J) 30,35,30 30 IA=I+(J*J-J)/2 ANORM=ANORM+A(IA)*A(IA) 35 CONTINUE IF(ANORM) 165,165,40 40 ANORM=1.414*DSQRT(ANORM) ANRMX=ANORM*RANGE/FLOAT(N) C C INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR C IND=0 THR=ANORM 45 THR=THR/FLOAT(N) 50 L=1 55 M=L+1 C C COMPUTE SIN AND COS C 60 MQ=(M*M-M)/2 LQ=(L*L-L)/2 LM=L+MQ 62 IF(DABS(A(LM))-THR) 130,65,65 65 IND=1 LL=L+LQ MM=M+MQ X=0.5*(A(LL)-A(MM)) 68 Y=-A(LM)/DSQRT(A(LM)*A(LM)+X*X) IF(X) 70,75,75 70 Y=-Y 75 SINX=Y/DSQRT(2.0*(1.0+(DSQRT(1.0-Y*Y)))) SINX2=SINX*SINX 78 COSX=DSQRT(1.0-SINX2) COSX2=COSX*COSX SINCS =SINX*COSX C C ROTATE L AND M COLUMNS C ILQ=N*(L-1) IMQ=N*(M-1) DO 125 I=1,N IQ=(I*I-I)/2 IF(I-L) 80,115,80 80 IF(I-M) 85,115,90 85 IM=I+MQ GO TO 95 90 IM=M+IQ 95 IF(I-L) 100,105,105 100 IL=I+LQ GO TO 110 105 IL=L+IQ 110 X=A(IL)*COSX-A(IM)*SINX A(IM)=A(IL)*SINX+A(IM)*COSX A(IL)=X 115 IF(MV-1) 120,125,120 120 ILR=ILQ+I IMR=IMQ+I X=R(ILR)*COSX-R(IMR)*SINX R(IMR)=R(ILR)*SINX+R(IMR)*COSX R(ILR)=X 125 CONTINUE X=2.0*A(LM)*SINCS Y=A(LL)*COSX2+A(MM)*SINX2-X X=A(LL)*SINX2+A(MM)*COSX2+X A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) A(LL)=Y A(MM)=X C C TESTS FOR COMPLETION C C TEST FOR M = LAST COLUMN C 130 IF(M-N) 135,140,135 135 M=M+1 GO TO 60 C C TEST FOR L = SECOND FROM LAST COLUMN C 140 IF(L-(N-1)) 145,150,145 145 L=L+1 GO TO 55 150 IF(IND-1) 160,155,160 155 IND=0 GO TO 50 C C COMPARE THRESHOLD WITH FINAL NORM C 160 IF(THR-ANRMX) 165,165,45 C C SORT EIGENVALUES AND EIGENVECTORS C 165 IQ=-N DO 185 I=1,N IQ=IQ+N LL=I+(I*I-I)/2 JQ=N*(I-2) DO 185 J=I,N JQ=JQ+N MM=J+(J*J-J)/2 IF(A(LL)-A(MM)) 170,185,185 170 X=A(LL) A(LL)=A(MM) A(MM)=X IF(MV-1) 175,185,175 175 DO 180 K=1,N ILR=IQ+K IMR=JQ+K X=R(ILR) R(ILR)=R(IMR) 180 R(IMR)=X 185 CONTINUE C RETURN END C_____________________________________________________________________________ C SUBROUTINE SORTP(NSTART,N) c c This routine is based on the SORT2 'heapsort' routine from Numerical c Recipes and sorts the quantities in vector WK from WK(NSTART) to WK(N) C in ascending order of magnitude and also accordingly rearranges vector C IPT of pointers to original positions of sorted quantities. c PARAMETER (NATMAX=500,NIBOND=4*NATMAX) c COMMON /sort/WK,ipt INTEGER*4 IPT(nibond),IIPT,L,I,J,IR REAL*8 WK(nibond),WWK C L=N/2+1 IR=N 10 CONTINUE IF(L.GT.NSTART)THEN L=L-1 WWK=WK(L) IIPT=IPT(L) ELSE WWK=WK(IR) IIPT=IPT(IR) WK(IR)=WK(1) IPT(IR)=IPT(1) IR=IR-1 IF(IR.EQ.NSTART)THEN WK(1)=WWK IPT(1)=IIPT RETURN ENDIF ENDIF I=L J=L+L 20 IF(J.LE.IR)THEN IF(J.LT.IR)THEN IF(WK(J).LT.WK(J+1))J=J+1 ENDIF IF(WWK.LT.WK(J))THEN WK(I)=WK(J) IPT(I)=IPT(J) I=J J=J+J ELSE J=IR+1 ENDIF GO TO 20 ENDIF WK(I)=WWK IPT(I)=IIPT GO TO 10 c RETURN END C C_____________________________________________________________________________ C subroutine cleans(line,leng,n) c c Correct to standard notation a numerical value written by FORTRAN c to first LENG characters of string LINE and return in N the position of c the first character c character line*80 c do 1 n=1,leng if(line(n:n).ne.' ')goto 2 1 continue c 2 if(line(n:n).eq.'.')then line(n-1:n-1)='0' n=n-1 return endif c if(line(n:n).eq.'-'.and.line(n+1:n+1).eq.'.')then line(n:n)='0' line(n-1:n-1)='-' n=n-1 return endif c return end C C_____________________________________________________________________________ c integer*2 function INKEY(N2) msps c c By L Pszczolkowski: c c...This emulates for MSF PS1.0 the INKEY function of Z.Czumaj which c in turn emulated for MSF5.0 the INKEY function from IIUWGRAF graphics c library for the Hercules card c c The function GETCHARQQ returns the ASCII character if the corresponding c key was pressed. If function or direction key was pressed then 0 c or hex E0 is returned and another call to GETCHARQQ is required to c get the extended code of the character c INCLUDE 'FLIB.FD' PS1 c USE MSFLIB PS4 c INTEGER*2 IK,n2 msps CHARACTER*1 KK msps c KK=GETCHARQQ() msps IK=ICHAR(KK) msps IF(IK.EQ.0 .OR. IK.EQ.224 ) THEN msps KK=GETCHARQQ() msps IK=-ICHAR(KK) msps ENDIF msps n2=ik msps INKEY=IK msps END msps C_____________________________________________________________________________ C_____________________________________________________________________________