C
C
C
C
C
C    NUPARM - A Nucleic Acid Parameter Determination Program 
C 
C                  M.Bansal & D.Bhattacharyya, 
C
C                  Molecular Biophysics Unit, 
C
C                 Indian Institute of Science, 
C
C                   Bangalore 560012, INDIA.
C
C
C   e-mail:   uunet!shakti!turing!mbu
C
C
C
C
C
C	     The NUPARM program has been developed so as to meet most of 
C
C	The requirements stipulated at the EMBO Workshop on DNA Curvature 
C
C	end  Bending,  held  at  Cambridge, U.K.,  in  Sept.  1988.   The 
C
C	Nomenclature  and  Description of the DNA  structural  parameters 
C
C	follows  the Cambridge Convention (as published in J. Mol.  Biol. 
C
C	(1989) vol. 205, 787-791, J. Biomol. Struct. Dynam. vol. 6,  627-
C
C	634(1989),  EMBO  J.  vol. 8, 1-4  (1989)).   While  considerable 
C
C	effort has been made to make NUPARM similar to the HELIB  program 
C
C	of  Dickerson et al, some differences are inevitable, because  of 
C
C       the  differences  in the procedure used.  We have also  made  the 
C
C	program  input  interactive rather than file  controlled.  For  a 
C
C	standard  PDB input file, no additional input file  is  required, 
C
C	and  the queries are only with regard to the choice of  base-pair 
C
C	center, Y-axis and global helix axis. 
C
C	
C     	      The NUPARM program creates three\four output files for each 
C	
C	run using an input file named "XXXX.YYY": 
C
C	1. A file "RUN.ANS" contains all the questions asked, as well  as 
C
C	the replies given, for later reference. 
C
C	2.  "XXXX.PRM" contains all the parameters calculated. 
C
C	3.  "XYZAXES.OUT" contains some additional information about  the 
C
C	local X, Y and Z axes at each step and the base-pair centers. 
C
C	4.  A  fourth  file  "XXXX.COOR" is  created  if  the  option  of 
C
C	reorienting the molecule about the global helix axis is used.  It 
C
C	contains the transformed coordinates, as well as related  residue 
C
C	and atom information. It can also be used as input for  NUPARM.   
C
C	
C   	      The definition of the local helix and wedge parameters  are 
C
C	in terms of the local helix axis and the mean Z-axis respectively 
C
C	for  the  doublet  involved (as described by  us  in  J.  Biomol. 
C
C	Struct. Dynam. (1989), vol. 6, pp.635-653). 
C
C	
C    	      There are several options available for obtaining the  best 
C
C	global  axis, but we have deliberately not allowed the option  of  
C
C	choosing  a few of these reference points and neglecting  others, 
C
C	as  in  HELIB/NEWHELIX.  Since the points chosen are  not  always 
C
C	clearly  defined  when  the structure  is  reported,  it  becomes 
C
C	difficult  for  other  workers  to make  a  clear  and  objective 
C
C	comparison between different structures. 
C
C	
C  	      Only  the  purine\pyrimidine ring atoms are  included  when 
C	
C	calculating the base normals.  We found only a nominal difference 
C
C	in  the parameters when all the pendant atoms are included.   The 
C
C	mean  base-pair  normal  is taken to be the average  of  the  two 
C
C	normals  in the base-pair, in order to minimize  the  differences 
C
C	due to the size of the purine and pyrimidine rings. 
C
C	
C             NUPARM  runs  interactively and asks  questions  which  are 
C	
C	meant to be self explanatory (we hope) and requires only a 'Y' or 
C
C	'N'  typed  on  the keyboard followed  by  <CR>\<ENTER>.  A  <CR> 
C
C	without  'Y' or 'N' is taken to indicate a 'No'.  It  first  asks 
C
C	for  the file name containing the atomic coordinates and  whether 
C
C	the  coordinates are in the Brookhaven Protein Data  Bank  format 
C
C	(PDB).  If not, it asks for the actual file format and unit  cell 
C
C	parameters.  It  can thus handle any other type  of  data  files, 
C
C	including   the  CORELS  and  KONNERT  formats   or    fractional 
C
C	coordinates.  The  format  should be given  within  brackets  and 
C
C	specify  (i)  an  atom name identifier (in  A4  format),  (ii)  a 
C
C	residue  (base) name identifier (in A3 format), (iii)  a  residue 
C
C	sequence  number identifier (in Integer format) and  (iv)  atomic 
C
C	coordinates.   The T-format specification should be used  if  the 
C
C	data for each atom does not follow the PDB order, viz. atom name, 
C
C	residue  name,  residue number and atomic coordinates.  The  atom 
C
C	name identifiers in the nucleotide residue should follow the  PDB 
C
C	code,  but  the  sequence  of atoms within  the  residue  is  not 
C
C	important. 
C
C
C	      The  program  next  asks whether  the  molecule  is  double 
C
C	stranded.  It  then  asks whether the  nucleotides  are  numbered 
C
C	according  to the standard PDB 5'-->3' residue numbering  scheme, 
C
C	with first residue paired to the last residue ( in the case of  a 
C
C	duplex).  If not, (eg. as in AMBER output, partially paired  tRNA 
C
C	or  any  other  type of structure) the file  containing  a  user-
C
C	defined residue numbering scheme has to be specified. The residue 
C
C	numbers  of  paired bases have to be given in 2I3  format  for  a 
C
C	duplex. 
C
C	
C   	      The  5'-->3'  direction of strand 1 is used to  define  the 
C	
C	positive Z-direction and the Y-axis is taken as pointing  towards 
C
C	this strand.  The X-axis then points towards the major groove, as 
C
C	suggested in the Cambridge Convention.  The base-pair center  can 
C
C	be  the center of gravity of the base-pair or the midpoint of  C6 
C
C	and  C8 atoms of pyrimidines and purines respectively.  The  long 
C
C	axis  (Y-axis)  can  be  along  the  C6---C8  or  the   C1'---C1' 
C
C	directions, but it always passes through the base-pair center. 
C
C	
C   	      For  a  double stranded structure, it asks  whether  single 
C
C	stranded  parameters  are  also required.  If  yes,  the  program 
C
C	calculates single strand parameters using a base as the repeating 
C
C	unit  (Local, local Helical and Global parameters) for  both  the 
C
C	strands separately. 
C
C	 
C    	      The   helix  can  be  reoriented  to  give  the   base-pair 
C
C	parameters  with respect to a global axis, as in  HELIB\NEWHELIX. 
C
C	The  realignment  is done in two steps:  (1)  The  best  line  is 
C
C	fitted   (by   the  method  of  least  squares)  to  a   set   of 
C
C	representative points taken either from all the residues, or  the 
C
C	doublet  steps.  The  molecule is rotated  such  that  this  line 
C
C	becomes   parallel   to  the  global  Z-axis.    (2)   The   same 
C
C	representative  points  are used to find the mean  point  through 
C
C	which  the line passes, and the molecule is translated such  that 
C
C	this  point coincides with the origin and the best line with  the 
C
C	Z-axis.   There  are several options available for  choosing  the 
C
C	representative points.  The user may choose between (i) the local 
C
C	helix  origins for each of the doublet steps, (ii) the  base-pair 
C
C	centers or (iii) any atom in the backbone.  
C
C	
C    	      The parameters calculated are grouped into four sets : 
C
C	1)  The  local base-pair step parameters :  Tilt(Tau),  Roll(Ro), 
C
C	Twist(Omega), Shift(Dx), Slide(Dy) and Rise(Dz). 
C
C	2)  The  local helical parameters (w.r.t. a helix axis  for  each 
C
C	doublet    step)   :   Inclination(Eta),   Tip(Theta),    Helical 
C
C	Twist(Omega), Displacements (dx and dy) and Helical Rise(dz). 
C
C	3)  The global helical parameters (w.r.t. the global  helix  axis 
C
C	for  the  whole molecule): Inclination(Eta),  Tip(Theta),  global 
C
C	Helical  Twist(Omega),  Displacements  (dx  and  dy)  and  global 
C
C	Helical Rise (dz). 
C
C	4)   Intra   base-pair  parameters  :   Propeller   Twist   (Pi), 
C
C	Buckle(Kappa),  Opening angle (defined as the angle  between  the 
C
C	lines  C1'---C8 and C1'---C6 projected on to the  mean  base-pair 
C
C	plane),  as well as the angles between the two  glycosidic  bonds 
C
C	(N---C1') and  the C1'---C1' vector. The  C8---C6  and  C1'---C1' 
C
C	distances   are  also  listed.  All  these  parameters  are   not 
C
C	calculated for a single strand molecule. 
C
C	
C   	      The direction cosines of the local helix axes and base-pair 
C	
C	normals  as  well as the coordinates of local helix  origins  and 
C
C	base-pair centers are  printed next. 
C
C	
C   	      The angle between each local helix axis and all others,  as 
C	
C	also the angles between base-pair normals are then tabulated. 
C
C	
C              An option is  available to obtain  the  cylindrical  polar 
C	
C	coordinates  of the phosphate atoms, assuming the global axis  to 
C	
C	be the Z-axis.  The intra and interchain phosphate distances  are 
C	
C	also listed.        
C
C
C	      An  option  is also available to calculate the  main  chain 
C
C	torsion  angles,  the glycosidic torsion angles and the  furanose 
C
C	ring,  endocyclic torsion angles.  The sugar pucker  phase  angle 
C
C	and   amplitude  of  puckering  (as  described  in   Harvey   and 
C
C	Prabhakaran, J. Amer. Chem. Soc. (1986), vol. 108,  pp.6128-6136) 
C
C	are listed  along with the ring torsion angles.  This calculation 
C
C	is the most time consuming part of the analysis program. 
C
C
C    
C    If any bugs/errors are detected please inform the authors.
C
C-----------------------------------------------------------------------
C
      DIMENSION IBPAIR(200),IBAPR(200),IPAIR(4),IP(4),IN(4),DC(3),
     1           DCR(3),RMAT(3,3),TMP(3),TITLE1(20)
	COMMON /CENTER/CENTRS(99,3)
      COMMON/ATINF/CR(2400,3),SECSQ(300),
     1             ISTS(300),IENS(300),NSEG,NAT
	COMMON /ATINFC/ATNM(2400),SECNM(300)
      COMMON/ATOMS/IBASE1(10),IBASE2(10),IBASE3(10),IBASE4(10)
      COMMON/GLOB/TILTG(99),ROLLG(99),TWISTG(99),DXG(99),DYG(99),HG(99)
      COMMON/LOCAL/TILTL(99),ROLLL(99),TWISTL(99),SLY(99),SLX(99),
     1DZL(99),LNDEX,DIETLT(99),DIEROL(99),PROTW(99),BUCAN(99),SLZ(99),
     2TWISTH(99),OPENAN(99),OPENDS(99),ANGBN9(99),OPENC1(99),ANGBN1(99)
      DOUBLE PRECISION XAXIS,YAXIS,ZAXIS,XM,YM,ZM1,ZM2,ZM3,Y1,Y2,Y3,
     1   Y12,Y22,Y32,X11,X12,X13,BMN,TL2,RL2,Y11,Y13
      COMMON /DISPL/DXLOC(99),DYLOC(99),DZLOC(99),KKTI,
     1      BSCENTR(99,3),HORIG(99,3)
	COMMON /BASINF/BASE(99)
      DOUBLE PRECISION ZAXISH,DC,DCR1,DCR2,DCR3,ANGL,ZGX,ZGY,ZGZ
      COMMON /ZSTRAXS/ZAXISH(3,99)
      COMMON /GLOREO/DCRS1,DCRS2,DCRS3,ANGLE,DX1,DY1,INFBP(2,99)
      COMMON /REPLY/ANSOR,ANSC1,ANSBPN,ANSDB,ANS,ANSHO
       CHARACTER *80 FILENM
       CHARACTER *40 FORMA,INPUT
       CHARACTER*80 TITLE1,TITLE
       CHARACTER*4  ATOM,ANORIN,ATNM
       CHARACTER*3 BAS1,BAS2,BAS3,BAS4,BASE,SECNM
       CHARACTER *1 ANS,ANSOR,ANSC1,ANSHO,ANSBPN,ANSPDB,ANSPP,ANSTOR,
     1ANSDB,ANSSNG
       DATA MAXSEG,MAXAT/99,1500/
C --  This programme can handle a maximum of 99 residues and 2400 atoms.
C
	ANORIN = '    '
      CONV=180.0/3.141592654
         INOHP = 0
         NB = 1
         IND = 0
C
	WRITE(6,9990)
 9990	FORMAT(1X,40('**'),//,5X,'NUPARM - A Nucleic Acid Parameter Dete
     1mination Program'//15X,'M. Bansal and D. Bhattacharyya',/18x,'Mole
     2cular Biophysics Unit',/17X,'Indian Institute of Science',/18X,'Ba
     3ngalore 560012,  INDIA',//5X,'e-mail:  uunet!shakti!turing!mbu')
C 
       OPEN(UNIT=8,FILE='RUN.ANS',STATUS='NEW')
       OPEN(UNIT=9,FILE='XYZAXES.OUT',STATUS='NEW')       
C
         WRITE(6,3)
  3     FORMAT(/' All the questions & answers are written into a file
     1 "RUN.ANS" ')
131      WRITE(6,4)
         WRITE(8,4) 
  4     FORMAT(/' Type in the name of the INPUT DATA file:')
	 READ(5,5) FILENM
  5     FORMAT(A80)
         WRITE(8,'(2X,''COORDINATES TAKEN FROM :''/1X,A80)') FILENM 
        OPEN(UNIT=4,FILE=FILENM,STATUS='OLD',ERR=131)
C
        CALL GETCOR(4,IERR,ANSPDB)
C
        IF(ANSPDB.EQ.'Y'.OR.ANSPDB.EQ.'y') THEN
	  REWIND 4
          IT = 0
          KKTI = 0
          DO WHILE(IT.EQ.0)
             READ(4,'(A80)') TITLE 
             IF(TITLE(1:6).EQ.'HEADER') THEN
	          KKTI = KKTI + 1
	          TITLE1(KKTI) = TITLE
             ELSE IF(TITLE(1:6).EQ.'COMPND') THEN
                 KKTI = KKTI + 1
                 TITLE1(KKTI) = TITLE
             ELSE IF(TITLE(1:6).EQ.'AUTHOR') THEN
                 KKTI = KKTI + 1
                 TITLE1(KKTI) = TITLE
	     ELSE IF(TITLE(1:4).EQ.'JRNL') THEN
	         KKTI = KKTI + 1
	         TITLE1(KKTI) = TITLE
             END IF
             IF(TITLE(1:4).EQ.'ATOM') IT = 1
             IF(KKTI.GT.19) IT = 1
          END DO
      END IF
C
      WRITE(6,91)
 91   FORMAT(/' Is the molecule a double stranded helix? (Y/N)')
      WRITE(8,91)
      READ(5,'(A1)')ANSDB
      IF(ANSDB.EQ.'y') ANSDB = 'Y'
      WRITE(8,'(5X,A1)')ANSDB
      NDB = NSEG/2
C
      WRITE(6,92)
 92   FORMAT(/' Are the residues numbered according to the Brookhaven fo
     1rmat? (Y/N)')
      WRITE(8,92)
      READ(5,'(A1)') ANSBPN
      IF(ANSBPN.EQ.'y') ANSBPN = 'Y'
      WRITE(8,'(5X,A1)') ANSBPN
      IF(ANSBPN .EQ. 'Y'.AND.ANSDB.EQ.'Y') THEN
           KOUNT = -1
	   DO KKK=1,NDB
              KOUNT = KOUNT + 2
              IBAPR(KOUNT) = KKK
	      IBAPR(KOUNT + 1) = NSEG - KKK + 1
	   END DO
	   NB = NSEG
       END IF
C
       IF(ANSBPN .NE. 'Y') THEN

133        WRITE(6,93)
 93   FORMAT(/' Give the file-name containing base-pair numbering inform
     1ation (in 2I3 Format):')
	  WRITE(8,93)
          READ(5,'(A40)') INPUT
       WRITE(8,'(2X,''Base-pair residue numbers from:''/1X,A40)')INPUT
	  OPEN(UNIT=2,FILE=INPUT,STATUS='OLD',ERR=133)
          DO  KB = 1,MAXSEG
                READ(2,10,END=1001) IB1,IB2
  10       FORMAT(2I3)
                IBAPR(NB) = IB1
                IBAPR(NB + 1) = IB2
                NB = NB + 2
          END DO
 1001    CONTINUE
          NB = NB - 1
        END IF
C
        IF(ANSDB .NE. 'Y' .AND.ANSBPN.EQ.'Y') THEN
            KOUNT = 1
           DO KK = 1,NSEG
              IBAPR(KOUNT) = KK
              IBAPR(KOUNT + 1) = 0
              KOUNT = KOUNT + 2
           END DO
          NB = NSEG * 2 
        END IF
C
          IF(IERR.EQ.1) THEN
             WRITE(6,*) 'ERROR IN DATA FILE'
             STOP
          END IF
C-----------------------------------------------------------------------
        CLOSE(UNIT=4)
C-----------------------------------------------------------------------
        ANSOR ='N'
C-- OPTION TO USE 'Centre of Gravity' AS ORIGIN IS SUPPRESSED.  
C   CAN BE ACTIVATED, "AT YOUR PERIL", BY DE-COMMENTING NEXT SIX LINES.
	IF(ANSDB.EQ.'Y') THEN
           WRITE(6,96)
 96   FORMAT(/' Should the base-pair ORIGIN be at the Centre of Gravity?  
     1'/' (NOT RECOMMENDED!! Default is midpoint of C6--C8:)(Y/N)')
           WRITE(8,96)
           READ(5,'(A1)') ANSOR
           IF(ANSOR.EQ.'y') ANSOR = 'Y'
	   WRITE(8,'(5X,A1)') ANSOR
	END IF
C-----------------------------------------------------------------------
C	IF THE MOLECULE IS SINGLE STRANDED THIS QUESTION IS OF NO USE

	IF ((ANSDB .EQ. 'Y')) THEN
         WRITE(6,97)
 97   FORMAT(/' Should the Y-axis be along the line joining C1'' atoms?' 
     1 /' (If "No",Y-axis is taken along the C6--C8 direction:)(Y/N)')
	 READ(5,'(A1)') ANSC1
         IF(ANSC1.EQ.'y') ANSC1 = 'Y'
	ELSE
	 ANSC1 = 'N'
	ENDIF

	 WRITE(8,97)
	 WRITE(8,'(5X,A1)') ANSC1
C
C-----------------------------------------------------------------------
        NRUN =1
	IF(ANSDB.EQ.'Y') THEN
         WRITE(6,'(/'' Are Single Strand Parameters required? (Y/N)'')')
	 READ(5,'(A1)') ANSSNG
         WRITE(8,'(/'' Are Single Strand Parameters required? (Y/N)'')')
	 WRITE(8,'(5X,A1)') ANSSNG
 	 IF(ANSSNG.EQ.'y'.OR.ANSSNG.EQ.'Y') NRUN = 3
	END IF
C
       DO 700  NTIM = 1,NRUN
          NG = 1
       IF(NTIM .EQ. 1)NG = 3
C
       DO  NGLOB = 1,NG
           LNDEX = 0
           NSTEPS = 0      
         DO 600 IBP = 1,(NB-2),2
             NSTEPS = NSTEPS + 1
            IF(IBP.EQ.NB-1) INOHP = 1
            IBAPR1 = IBAPR(IBP)
            IBAPR2 = IBAPR(IBP + 1)
          IF(NTIM .EQ. 2)IBAPR2 = 0
            IF(NTIM .EQ. 3) THEN
               IBAPR1 = IBAPR(NB - IBP + 1)
               IBAPR2 = 0
            END IF
C
         CALL PICATP(IBAPR1,IBAPR2,IPAIR,IBASE1,IBASE2,NB1,NB2,INOHP,
     1           BAS1,BAS2)

            BASE(NSTEPS)(1:1) = BAS1(1:1)
            BASE(NSTEPS)(2:2) = ':'
            BASE(NSTEPS)(3:3) = BAS2(1:1)
	    INFBP(1,NSTEPS) = IBAPR1
	    INFBP(2,NSTEPS) = IBAPR2
C 
            IP(1) = IPAIR(2)
            IP(2) =  IPAIR(4)
            IN(1) = IPAIR(1)
            IN(2) = IPAIR(3)
C
            IBAPR3 = IBAPR(IBP+2)
            IBAPR4 = IBAPR(IBP + 3)
             IF(NTIM .EQ. 2)IBAPR4 = 0
               IF(NTIM .EQ. 3) THEN
                  IBAPR3 = IBAPR(NB - IBP -1)
                  IBAPR4 = 0
               END IF
C
        CALL PICATP(IBAPR3,IBAPR4,IPAIR,IBASE3,IBASE4,NB3,NB4,INOHP,
     1           BAS3,BAS4)
            BASE(NSTEPS+1)(1:1) = BAS3(1:1)
            BASE(NSTEPS+1)(2:2) = ':'
            BASE(NSTEPS+1)(3:3) = BAS4(1:1)
	    INFBP(1,NSTEPS+1) = IBAPR3
	    INFBP(2,NSTEPS+1) = IBAPR4
            IP(3) = IPAIR(2)
            IP(4) = IPAIR(4)
            IN(3) = IPAIR(1)
            IN(4) = IPAIR(3)
       CALL FINDVEC(IP,IN,NB1,NB2,NB3,NB4,NSTEPS,ANSOR,ANSC1,ANS,NTIM)
 600  CONTINUE
C
      IF(NTIM .GT. 1)GO TO 500
        IF(NGLOB .EQ. 1) THEN
         WRITE(6,99)
         WRITE(8,99)
 99   FORMAT(/,' Should the molecule be reoriented about a LINEAR GLOBAL
     1 Helix Axis? (Y/N)')
        READ(5,'(A1)')ANS
	IF(ANS.EQ.'y') ANS = 'Y'
        WRITE(8,'(5X,A1)')ANS
        IF(ANS .NE.'Y') GO TO 500
C     
        END IF
C
           NPTS = NSTEPS
         IF(NGLOB .EQ. 1)THEN
	REWIND 9
	   CALL REORNT(IBAPR,NSTEPS,ANORIN,NPOINT)
	   IF(ANORIN.EQ.'O   '.OR.ANORIN.EQ.'C   ') THEN

 	      CALL LINEFIT(DC,NPOINT)
	   ELSE IF(ANSDB.EQ.'Y') THEN
	      CALL LINETM(DC,NPOINT)
	   ELSE
	      CALL LINESN(DC,NPOINT)
	   END IF
	   ZGX = 0.0000000000
	   ZGY = 0.0000000000
   	   ZGZ = 1.0000000000
           CALL CROSS(DC(1),DC(2),DC(3),ZGX,ZGY,ZGZ,DCR1,DCR2,DCR3)
           ANGL = DACOS(DC(3))
           ANGLE = ANGL * CONV
           DCRS1 = DCR1 * 1.0
           DCRS2 = DCR2 * 1.0
           DCRS3 = DCR3 * 1.0
           CALL ROTMAT(DCRS1,DCRS2,DCRS3,ANGLE,RMAT)
            DO KJ=1,NAT
              CALL MATMUL(RMAT,CR(KJ,1),CR(KJ,2),CR(KJ,3),TMP)
                DO KL=1,3
                     CR(KJ,KL) = TMP(KL)
                END DO
            END DO
C
          ELSE IF(NGLOB .EQ. 2)THEN    
	         REWIND 9
 	         CALL REORNT(IBAPR,NPTS,ANORIN,NPOINT)
	   IF(ANORIN.EQ.'O   '.OR.ANORIN.EQ.'C   ') THEN
                 CALL  POINT(NPOINT,DX1,DY1)
	   ELSE IF(ANSDB.EQ.'Y') THEN
	      CALL POINDB(NPOINT,DX1,DY1)
	   ELSE
	      CALL POINSN(NPOINT,DX1,DY1)
	   END IF
                 DO KJ=1,NAT
                     CR(KJ,1) = CR(KJ,1) - DX1
                     CR(KJ,2) = CR(KJ,2) - DY1
                 END DO
           END IF
C 
        END DO     
 500  CONTINUE
        CALL WRITING(FILENM,TITLE1,INPUT,NTIM,NRUN,ANORIN)
C-----------------------------------------------------------------------
      IF(NTIM .GT. 1)GO TO 701
      IF(ANS .EQ. 'Y')CALL REWRIT(FILENM,TITLE1,KKTI)
C   Reoriented atomic coordinates written in xxx.COOR
C-----------------------------------------------------------------------
	IF(NTIM.EQ.1.AND.ANS.EQ.'Y') WRITE(9,1101)
 701	CONTINUE
	IF(NTIM.EQ.2.AND.ANS.EQ.'Y') WRITE(9,1102)
 1101	FORMAT('1',30X,'First strand axes etc.'/20X,3('--------------'))
 1102	FORMAT('1',30X,'Second strand axes etc.'/20X,3('-------------'))
 700  CONTINUE
C-----------------------------------------------------------------------
      WRITE(6,9993)
      WRITE(8,9993)
9993  FORMAT(/' Are P--P distances and Cylindrical Polar Coordinates
     1 required?(Y/N)' /' Only meaningful for reoriented molecule.')
      READ(5,'(A1)')ANSPP
      IF(ANSPP.EQ.'y') ANSPP = 'Y'
      WRITE(8,'(5X,A1)')ANSPP
      IF(ANSPP.EQ.'Y' )CALL RTHEPH(ANSDB)
C-----------------------------------------------------------------------
       WRITE(6,9995)
       WRITE(8,9995)
9995  FORMAT(/' Is the torsion angle calculation required? (Y/N)')
       READ(5,'(A1)')ANSTOR
       IF(ANSTOR.EQ.'y') ANSTOR = 'Y'
       WRITE(8,'(5X,A1)')ANSTOR
       IF(ANSTOR.EQ.'Y') CALL TORSION(ANSDB,ANSBPN)
C-----------------------------------------------------------------------     
       STOP
       END
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------  
      SUBROUTINE GETCOR(IUNIT,IERR,ANSPDB)
C..
C     THIS ROUTINE IS TAKEN FROM THE PROGRAM CHEM
C     AND MODIFIED
C..   READING COORDINATES FROM IUNIT IN PDB FORMAT
C
      COMMON/ATINF/CR(2400,3),SECSQ(300),ISTS(300)
     1            ,IENS(300),NSEG,NAT
	COMMON /ATINFC/ATNM(2400),SECNM(300)
 
C
      CHARACTER*1 ANSPDB
      CHARACTER*80 LINE
      CHARACTER*3 RESN,SECNM
      CHARACTER*4 ATNM
      CHARACTER*40 FORMA
      INTEGER RESNO
      DATA MAXSEG,MAXAT/300,2400/
 5    FORMAT(A80)
 15   FORMAT(12X,A4,1X,A3,2X,I4,4X,3F8.3)
C
  6   FORMAT('IO ERROR READING *PDB* FILE')
 16   FORMAT('SYSTEM TOO LARGE')
 26   FORMAT('BAD ORDERING OF RESIDUES')
 36   FORMAT('NO ATOMS DEFINED WITHIN *PDB* FILE')
C
      NW = 6
      REWIND IUNIT
C
      IERR = 0
      NAT = 0
      NSEG = 0
      LRESNO = -1
	WRITE(6,996)
	WRITE(8,996)
 996  FORMAT(/' Is the data file in the "BROOKHAVEN PDB" format?(Y/N)')
	READ(5,'(A1)') ANSPDB
        IF(ANSPDB.EQ.'y') ANSPDB = 'Y'
	WRITE(8,'(5X,A1)') ANSPDB
	IF(ANSPDB.EQ.'Y') THEN
C
C     LOOP TO READ THE COORDINATES IN PDB FORMAT
C
 10   CONTINUE
          READ(4,5,END=61,ERR=21) LINE
          IF (LINE(1:4).EQ.'ATOM') THEN
              NAT = NAT + 1
              IF (NAT.GT.MAXAT) GO TO 31
              READ(LINE,15) ATNM(NAT),RESN,RESNO,(CR(NAT,L),L=1,3)
              IF(ATNM(NAT)(1:1).EQ.' ') THEN
                 ATNM(NAT)(1:1)=ATNM(NAT)(2:2)
                 ATNM(NAT)(2:2)=ATNM(NAT)(3:3)
                 ATNM(NAT)(3:3)=ATNM(NAT)(4:4)
                 ATNM(NAT)(4:4)=' '
              ENDIF
              IF (RESNO.NE.LRESNO) THEN
C                 DO 20 KSEG = 1,NSEG
C                 IF (SECSQ(KSEG).EQ.RESNO) GO TO 41
C20               CONTINUE
                  NSEG = NSEG + 1
                  IF (NSEG.GT.MAXSEG) GO TO 31
                  SECSQ(NSEG) = NSEG
                  SECNM(NSEG) = RESN
                  ISTS(NSEG) = NAT
                  LRESNO = RESNO
              ENDIF
              IENS(NSEG) = NAT
          ENDIF
      GO TO 10
C     END OF THE READING LOOP
 61   CONTINUE
	ELSE
	    WRITE(6,997)
	    WRITE(8,997)
 997  FORMAT(/' If not, type in the input data format'/' Field specific  
     1ations required: Atom Name, Res.Name, Res.No., XYZ Coords.')
	    READ(5,'(A40)') FORMA
	    WRITE(8,'(A40)') FORMA
	    WRITE(6,998)
	    WRITE(8,998)
 998  FORMAT(/' TYPE UNIT CELL PARAMETERS: A,B,C AND ALPHA,BETA,GAMA')
C
	    READ(5,*) AC,BC,CC,ALP,BET,GAM
	    WRITE(8,911) AC,BC,CC,ALP,BET,GAM
 911  FORMAT(/2X,'UNIT CELL PARAMETERS :: A =',F8.2,' B =',F8.2,' C =',
     1       F8.2,/21X,' APLHA =',F8.2,' BETA =',F8.2,' GAMA =',F8.2)
 101	CONTINUE
              READ(4,5,END=161,ERR=21) LINE
C
              NAT = NAT + 1
              IF (NAT.GT.MAXAT) GO TO 31
              READ(LINE,FORMA) ATNM(NAT),RESN,RESNO,(CR(NAT,L),L=1,3)
              IF(ATNM(NAT)(1:1).EQ.' ') THEN
                 ATNM(NAT)(1:1)=ATNM(NAT)(2:2)
                 ATNM(NAT)(2:2)=ATNM(NAT)(3:3)
                 ATNM(NAT)(3:3)=ATNM(NAT)(4:4)
                 ATNM(NAT)(4:4)=' '
              ENDIF
C
              IF (RESNO.NE.LRESNO) THEN
                  NSEG = NSEG + 1
                  IF (NSEG.GT.MAXSEG) GO TO 31
                  SECSQ(NSEG) = NSEG
                  SECNM(NSEG) = RESN
                  ISTS(NSEG) = NAT
                  LRESNO = RESNO
              ENDIF
C
              IENS(NSEG) = NAT
      GO TO 101
C     END OF THE READING LOOP
 161  CONTINUE
	CALL FRACTION(AC,BC,CC,ALP,BET,GAM)
	END IF
      IF (NAT.EQ.0) GO TO 51
      IERR = 0
C      WRITE(6,*) 'Total no. of ATOMS in the "PDB" file ::',NAT
      RETURN
C
 21   CONTINUE
      WRITE (NW,6)
      IERR = 1
      RETURN
C
 31   CONTINUE
      WRITE (NW,16)
      IERR = 1
      RETURN
C
C 41   CONTINUE
C     WRITE (NW,26)
C
 51   CONTINUE
      WRITE (NW,36)
      IERR = 1
      RETURN
      END
 
      SUBROUTINE PICATP(IBAPR1,IBAPR2,IPAIR,IBASE1,IBASE2,NB1,
     1                  NB2,INOHP,BASE1,BASE2)
C
         COMMON/ATINF/CR(2400,3),SECSQ(300),
     1             ISTS(300),IENS(300),NSEG,NAT
	COMMON /ATINFC/ATNM(2400),SECNM(300)
 
C
      DIMENSION IPAIR(4),ATOMPA(9)
      DIMENSION IBASE1(10),IBASE2(10)
C
      CHARACTER*1 PRIME
      CHARACTER*3 SECNM,TMP
      CHARACTER*3 BASE1,BASE2
      CHARACTER*4 ATNM,ATOMPA
C
      DATA ATOMPA/'C4'' ','C3'' ','O4'' ','C2'' ','C1'' ','N9  ','N1  ',
     1'C8  ','C6  '/
      DATA PRIME/''''/
C   
      DO KL=1,4
      IPAIR(KL) = 0
      END DO
C
C     TO CHANGE THE C1* etc to C1' AND O1' TO O4' TO MATCH
C     THE ATOM NAMES FROM PDB AND AMBER
C
      DO 345 IA = 1, NAT
      IF(ATNM(IA)(1:3).EQ.'O1'' ')ATNM(IA)=ATOMPA(3)
      IF(ATNM(IA)(3:3).EQ.'*')ATNM(IA)(3:3)=PRIME
 345  CONTINUE
C     WRITE(6,555) NSEG,(SECNM(IR),IR=1,NSEG)
 555  FORMAT(3X,I5,30(2X,A1))

      DO 445 IR = 1,NSEG
      IF(SECNM(IR)(3:3).EQ.'G')SECNM(IR)='GUA'
      IF((SECNM(IR)(3:3).EQ.'A') .AND.
     1 (SECNM(IR) .NE. 'URA') .AND.
     1 (SECNM(IR) .NE.'GUA')) SECNM(IR)='ADE'
      IF(SECNM(IR)(3:3).EQ.'C')SECNM(IR)='CYT'
      IF((SECNM(IR)(3:3).EQ.'T') .AND.
     1 (SECNM(IR) .NE. 'CYT'))SECNM(IR)='THY'
      IF(SECNM(IR)(3:3).EQ.'U')SECNM(IR)='URA'
 445  CONTINUE

      ITYPE1 = 2
      ITYPE2 = 2
      IF(SECNM(IBAPR1).EQ.'ADE'.OR.SECNM(IBAPR1).EQ.'GUA') ITYPE1=1
      IF(SECNM(IBAPR2).EQ.'ADE'.OR.SECNM(IBAPR2).EQ.'GUA') ITYPE2=1
      ISTR = ISTS(IBAPR1)
      BASE1 = SECNM(IBAPR1)
      IEND = IENS(IBAPR1)
      DO 100 IAT = ISTR,IEND
      IF(ITYPE1.EQ.1) THEN
      IF(ATNM(IAT)(1:2).EQ.ATOMPA(6)(1:2)) IPAIR(2)=IAT
      IF(ATNM(IAT)(1:2).EQ.ATOMPA(8)(1:2)) IPAIR(1)=IAT
      ELSE IF(ITYPE1.EQ.2) THEN
      IF(ATNM(IAT)(1:2).EQ.ATOMPA(7)(1:2)) IPAIR(2)=IAT
      IF(ATNM(IAT)(1:2).EQ.ATOMPA(9)(1:2)) IPAIR(1)=IAT
      ENDIF
100   CONTINUE
C
      CALL BASPIC(ISTR,IEND,ITYPE1,IBASE1,NB1)
C
      IF(IBAPR2.NE.0) THEN
C
      ISTR = ISTS(IBAPR2)
      IEND = IENS(IBAPR2)
      BASE2 = SECNM(IBAPR2)
      DO 101 IAT = ISTR,IEND
      IF(ITYPE2.EQ.1) THEN
      IF(ATNM(IAT)(1:2).EQ.ATOMPA(6)(1:2)) IPAIR(4)=IAT
      IF(ATNM(IAT)(1:2).EQ.ATOMPA(8)(1:2)) IPAIR(3)=IAT
      ELSE IF(ITYPE2.EQ.2) THEN
      IF(ATNM(IAT)(1:2).EQ.ATOMPA(7)(1:2)) IPAIR(4)=IAT
      IF(ATNM(IAT)(1:2).EQ.ATOMPA(9)(1:2)) IPAIR(3)=IAT
      ENDIF
101   CONTINUE
C    
      CALL BASPIC(ISTR,IEND,ITYPE2,IBASE2,NB2)
C
      ELSE
C
      IPAIR(3) = 0
      BASE2 = '   '
C
      DO IAT=ISTR,IEND
      IF(ATNM(IAT)(1:3).EQ.ATOMPA(5)(1:3)) IPAIR(4) = IAT
      END DO
      ENDIF
C
      RETURN
      END

C----------------------------------------------------------------------
      SUBROUTINE BASPIC(ISTR,IEND,ITYPE,IATPIC,NB)
      COMMON/ATINF/CR(2400,3),SECSQ(300),
     1             ISTS(300),IENS(300),NSEG,NAT
	COMMON /ATINFC/ATNM(2400),SECNM(300)
C
C
      DIMENSION ATPUR(10),ATPYR(7),IATPIC(10)
      CHARACTER*3 SECNM
      CHARACTER*4 ATNM,ATPUR,ATPYR
 
      DATA ATPUR/'C1'' ','N9  ','C8  ','N7  ','C5  ','C6  ','N1  ',
     1            'C2  ','N3  ','C4  '/
      DATA ATPYR/'C1'' ','N1  ','C6  ','C5  ','C4  ','N3  ','C2  '/
      DATA NPUR,NPYR/10,7/
C
C     WRITE(9,887) ITYPE
887   FORMAT(' ITYPE ' , I4/)
      GO TO (20,30),ITYPE
C
C     PURINE  BASE ATOMS
C
20    CONTINUE
      NB=NPUR
      DO 100 IA = ISTR,IEND
      DO 100 IP = 1,NPUR
      IF(ATNM(IA)(1:3).EQ.ATPUR(IP)(1:3)) IATPIC(IP)=IA
100   CONTINUE
C     WRITE(9,57)(IATPIC(IC),IC=1,NPUR)
57    FORMAT(20I4)
      RETURN
30    CONTINUE
      NB=NPYR
      DO 200 IA = ISTR,IEND
      DO 200 IP = 1,NPYR
      IF(ATNM(IA)(1:3).EQ.ATPYR(IP)(1:3)) IATPIC(IP)=IA
200   CONTINUE
C     WRITE(9,57)(IATPIC(IC),IC=1,NPYR)
      RETURN
      END
 
C----------------------------------------------------------------------!
      SUBROUTINE FINDVEC(IPS,INS,NAT1,NAT2,NAT3,NAT4,NOBASE,ANSOR,ANSC1,
     1ANS,NTIM)
C----------------------------------------------------------------------!
      DOUBLE PRECISION EL,EM,EN,EL2,EM2,EN2,XAX,XAXIS,YAXIS,DOTPR,BM,BMN
     1 ,WGLT,WGLL,WGTILT,WGROLL,PROJ,ELT,EMT,ENT,ELN,EMN,ENN,ZAXIS,
     2 BPDIR,BPNORML,B,C8C6,YC8C6
      DOUBLE PRECISION EL1,EM1,EN1,SM,WTILT,WROLL,WGLT1,WGRL1
      DOUBLE PRECISION BNORM,EK1,EK2,EK3,VEC,TEMPX,TEMPY 
      DIMENSION X(3,15),ATOM(15),NSP(15),ORIG(3),PROJ(3),XYZ1(3),XYZ2
     1(3),XYZ3(3),XYZ4(3),DVA(15),N(15),PROJ1(2),B(2,3),AX(3),IPS(4),
     2C(3),XAX(3),XAXIS(2,3),YAXIS(2,3),BMN(3),INS(4),NOATM(2),VEC(3)
      DIMENSION X1(3,15),ATOM1(15),NSP1(15),N1(15),XC8(3),XC6(3),RC8(3)
     1         ,RC6(3),SM(3),BM(3),BNORM(3,4),ORIGIN(2,3),TEMPX(3),
     2         TEMPY(3),R(3,3),XC1P(3),RC1P(3),C8C1(3),C6C1(3),XC1P2(3),
     3ZAXIS(2,3),BPDIR(3,4),ORIG1(3),ORIG2(3),XC81(3),XC82(3),XC1P1(3)
     4         ,C8C6(3),YC8C6(3),XMID(3)
      COMMON /GLOB/TILTG(99),ROLLG(99),TWISTG(99),DXG(99),DYG(99),HG(99)
      COMMON /LOCAL/TILTL(99),ROLLL(99),TWISTL(99),SLY(99),SLX(99),
     1DZL(99),LNDEX,DIETLT(99),DIEROL(99),PROTW(99),BUCAN(99),SLZ(99),
     2TWISTH(99),OPENAN(99),OPENDS(99),ANGBN9(99),OPENC1(99),ANGBN1(99)
       COMMON /BPZAXS/BPNORML(99,3)
       COMMON /DISPL/DXLOC(99),DYLOC(99),DZLOC(99),KKTI,
     1       BSCENTR(99,3),HORIG(99,3)
	COMMON /BASINF/BASE(99)
       COMMON/ATINF/CR(2400,3),SECSQ(300),ISTS(300),IENS(300),NSEG,NAT
	COMMON /ATINFC/ATNM(2400),SECNM(300)
       COMMON /ATOMS/IBASE1(10),IBASE2(10),IBASE3(10),IBASE4(10)
      LOGICAL DOUBL1,DOUBL2,STKINV
      CHARACTER*4 ATOM,ATOM1,ATNM
      CHARACTER*3 BASE,SECNM
      CHARACTER*1 TMPBAS,ANSOR,ANSC1,ANS
      CONV=180.0/3.141592654
      NRTCAL = 0
C	REWIND 9
      DOUBL1 = .TRUE.
      STKINV = .FALSE.
      DOUBL2 = .TRUE.
C
C	WRITE(9,'(3(''****************************'')')
        DO  J=1,NAT1
                    IF(IBASE1(J).EQ.INS(1)) NC8 = J
                    IF(IBASE1(J).EQ.INS(2)) NC6 = J
                    IF(IBASE1(J).EQ.IPS(1)) NN9 = J
                    IF(IBASE1(J).EQ.IPS(2)) NN1 = J
            KJ = IBASE1(J)
            ATOM(J) = ATNM(KJ)
            N(J) = 0
           DO KC =1,3
              X(KC,J) = CR(KJ,KC)
           END DO
C     
         IF(ATOM(J).EQ.'C1'' ') THEN
            DO KC=1,3
               XC1P(KC) = X(KC,J)
	       XC1P1(KC) = X(KC,J)
            END DO
	 ELSE
           N(J) = J
         ENDIF
C
       END DO
C
        DO KC = 1,3
        XC8(KC) = CR(INS(1),KC)
	XC81(KC) = CR(INS(1),KC)
        XC6(KC) = CR(INS(2),KC)
        XYZ1(KC) = CR(IPS(1),KC)
        XYZ4(KC) = CR(IPS(2),KC)
        XYZ2(KC) = CR(INS(1),KC)
        XYZ3(KC) = CR(INS(2),KC)
        END DO
C
      IF(INS(2).NE.0) THEN
           DO  I=1,NAT2
                    IF(IBASE2(I).EQ.INS(1)) NC8 = I
                    IF(IBASE2(I).EQ.INS(2)) NC6 = I
                    IF(IBASE2(I).EQ.IPS(1)) NN9 = I
                    IF(IBASE2(I).EQ.IPS(2)) NN1 = I
            KJ = IBASE2(I)
            ATOM1(I) = ATNM(KJ)
            N1(I) = 0
           DO KC=1,3
              X1(KC,I) = CR(KJ,KC)
           END DO
C
        IF(ATOM1(I).EQ.'C1'' ') THEN
           DO KC=1,3
                    RC1P(KC) = X1(KC,I)
           END DO
	ELSE
           N1(I) = I
        ENDIF
      END DO
C --------------------------------------------------------------------
C   ASSIGNED THE ATOM NOS. FOR FIXING Y-AXIS .
C   NOW FINDING THE COORDINATES OF BASE PAIR ORIGIN AND DIRECTION
C   RATIOS OF THE "Y"-AXIS.
C --------------------------------------------------------------------
         DO KJ=1,3
           C8C1(KJ) = XC1P(KJ) - XC8(KJ)
           C6C1(KJ) = RC1P(KJ) - XC6(KJ)
         END DO
       SUMANG = 0.0
       SUMDS1 = 0.0
       SUMDS2 = 0.0
       DO KJ=1,3
         SUMDS1 = C8C1(KJ) * C8C1(KJ) + SUMDS1
         SUMDS2 = C6C1(KJ) * C6C1(KJ) + SUMDS2
       END DO
       SUMDS1 = SQRT(SUMDS1)
       SUMDS2 = SQRT(SUMDS2)
       DO KJ=1,3
          C8C1(KJ) = C8C1(KJ)/SUMDS1
          C6C1(KJ) = C6C1(KJ)/SUMDS2
      	  BPDIR(KJ,1) = C8C1(KJ) * 1.0
	  BPDIR(KJ,2) = C6C1(KJ) * 1.0
       END DO
	XN9C1 = 0.0
	XN1C1 = 0.0
	C1C1DS = 0.0
	C8C6DS = 0.0
	DSN1C1 = 0.0
	DSN9C1 = 0.0
       DO KJ=1,3
       C8C6DS = C8C6DS + (XC8(KJ) - XC6(KJ)) * (XC8(KJ) - XC6(KJ))
	C1C1DS = C1C1DS + (XC1P(KJ) - RC1P(KJ))*(XC1P(KJ) - RC1P(KJ))
	XN9C1 = XN9C1 + (X(KJ,NN9)-XC1P(KJ))*(RC1P(KJ)-XC1P(KJ))
	XN1C1 = XN1C1 + (X1(KJ,NN1)-RC1P(KJ))*(XC1P(KJ)-RC1P(KJ))
	DSN9C1 = DSN9C1 + (X(KJ,NN9)-XC1P(KJ))*(X(KJ,NN9)-XC1P(KJ))
	DSN1C1 = DSN1C1 + (X1(KJ,NN1)-RC1P(KJ))*(X1(KJ,NN1)-RC1P(KJ))
       END DO
	DSN9C1 = SQRT(DSN9C1)
	DSN1C1 = SQRT(DSN1C1)
        OPENDS(LNDEX+1) = SQRT(C8C6DS)
	OPENC1(LNDEX+1) = SQRT(C1C1DS)
	ANGBN9(LNDEX+1) = ACOS(XN9C1/(DSN9C1*OPENC1(LNDEX+1)))*CONV
	ANGBN1(LNDEX+1) = ACOS(XN1C1/(DSN1C1*OPENC1(LNDEX+1)))*CONV
C
       DO KJ=1,3
        IF(ANSOR.EQ.'Y') THEN
           ORIG1(KJ)=0.0
           ORIG2(KJ)=0.0
       	   DO NN=1,NAT1	    
             IF(ATOM(NN) .NE. 'C1'' ')ORIG1(KJ) = ORIG1(KJ) + X(KJ,NN)
	   END DO
             ORIG1(KJ) = ORIG1(KJ)/(NAT1-1)
	   DO NN=1,NAT2
	     IF(ATOM1(NN).NE.'C1'' ')ORIG2(KJ) = ORIG2(KJ) + X1(KJ,NN) 
	   END DO
	     ORIG2(KJ) = ORIG2(KJ)/(NAT2-1)
           ORIG (KJ) = (ORIG1(KJ) + ORIG2(KJ)) * 0.5
         ELSE
           ORIG(KJ) = (XC8(KJ) + XC6(KJ)) * 0.5
	 END IF
        ORIGIN(1,KJ) = ORIG(KJ)
        PROJ(KJ) =  XC8(KJ) - XC6(KJ)
       END DO
C
 	 IF(ANSC1.EQ.'Y') THEN
	   DO KKK=1,3
	      PROJ(KKK) = XC1P(KKK) - RC1P(KKK)
	      C8C6(KKK) = XC8(KKK) - XC6(KKK)
	      XMID(KKK) = (XC1P(KKK) + RC1P(KKK)) * 0.5
	   END DO
	   CALL NORMAL(C8C6(1),C8C6(2),C8C6(3),YC8C6(1),YC8C6(2),YC8C6(3
     1))
	   CALL NORMAL(PROJ(1),PROJ(2),PROJ(3),YAXIS(1,1),YAXIS(1,2),
     1YAXIS(1,3))
	   CONSTN=(YAXIS(1,1)*(XMID(1)-XC8(1))+YAXIS(1,2)*(XMID(2)-XC8(2
     1))+YAXIS(1,3)*(XMID(3)-XC8(3)))/(YAXIS(1,1)*YC8C6(1)+YAXIS(1,2)*
     2 YC8C6(2)+YAXIS(1,3)*YC8C6(3))
	XYZMD1 = XC8(1) + YC8C6(1) * CONSTN
	XYZMD2 = XC8(2) + YC8C6(2) * CONSTN
	XYZMD3 = XC8(3) + YC8C6(3) * CONSTN
	ELSE
	   CALL NORMAL(PROJ(1),PROJ(2),PROJ(3),YAXIS(1,1),YAXIS(1,2),
     1 YAXIS(1,3))
	END IF
	IF((ANSOR.EQ.'Y').OR.(ANSC1.NE.'Y')) THEN
	   ORIGIN(1,1) = ORIG(1)
	   ORIGIN(1,2) = ORIG(2)
	   ORIGIN(1,3) = ORIG(3)
	ELSE
	   ORIGIN(1,1) = XYZMD1
	   ORIGIN(1,2) = XYZMD2
	   ORIGIN(1,3) = XYZMD3
	END IF
C    
	WRITE(9,'('' BASE-PAIR ORIGIN-1'',3F10.4)')(ORIGIN(1,KJ),KJ=1,3)
       CALL NORMAL(PROJ(1),PROJ(2),PROJ(3),YAXIS(1,1),YAXIS(1,2),
     1 YAXIS(1,3))
C
C----------------------------------------------------------------------
C FIND L,M,N, OF THE TWO BASE NORMALS & B.P. NORMAL IN FIRST base-pair
C----------------------------------------------------------------------
C
       CALL PLANE(X,NAT1,N,1,DVA,EL,EM,EN,DETR,DVS1)
       BNORM(1,1) = EL
       BNORM(2,1) = EM
       BNORM(3,1) = EN
C
       CALL PLANE(X1,NAT2,N1,1,DVA,EL1,EM1,EN1,DETR,DVS2)
       ANG = EL * EL1 + EM * EM1 + EN * EN1
       IF(ANG.LT.0.0) THEN
         EL1 = -EL1
         EM1 = -EM1
         EN1 = -EN1
       ENDIF
        BNORM(1,2) = EL1
        BNORM(2,2) = EM1
        BNORM(3,2) = EN1
C
C MEAN OF THE BASE-NORMALS TO GET B.P. NORMAL
        EL1 = EL + EL1
        EM1 = EM + EM1
        EN1 = EN + EN1
        CALL NORMAL(EL1,EM1,EN1,EL,EM,EN)
C 
      ELSE
C
         IF(NC8.EQ.0.AND.NC6.NE.0) THEN
              NC8 = NC6
              DO KL = 1,3
                    XC8(KL) = XC6(KL)
              END DO
         ENDIF
      CALL SINGAXS(IPS,N,NAT1,NC8,X,1,XAXIS,YAXIS,EL,EM,EN,IBASE1,S1)
C
       DOUBL1 = .FALSE.
C
      ENDIF
C---------------------------------------------------------------------- 
C  SECOND BASE/base-pair OF THE DOUBLET
C----------------------------------------------------------------------
       DO J=1,NAT3
                    IF(IBASE3(J).EQ.INS(3)) NC8 = J
                    IF(IBASE3(J).EQ.INS(4)) NC6 = J
                    IF(IBASE3(J).EQ.IPS(3)) NN9 = J
                    IF(IBASE3(J).EQ.IPS(4)) NN1 = J
           KJ = IBASE3(J)
           ATOM(J) = ATNM(KJ)
            N(J) = 0
           DO KC = 1,3
              X(KC,J) = CR(KJ,KC)
           END DO
C   
         IF(ATOM(J).EQ.'C1'' ') THEN
              DO KC=1,3
                  XC1P(KC) = X(KC,J)
	          XC1P2(KC) = X(KC,J)
              END DO
         ELSE
           N(J) = J
	 END IF
       END DO
C
        DO KC = 1,3
        RC8(KC) = CR(INS(3),KC)
	XC82(KC) = CR(INS(3),KC)
        RC6(KC) = CR(INS(4),KC)
        XYZ1(KC) = CR(IPS(3),KC)
        XYZ4(KC) = CR(IPS(4),KC)
        XYZ2(KC) = CR(INS(3),KC)
        XYZ3(KC) = CR(INS(4),KC)
        END DO
C
      IF(INS(4).NE.0) THEN
C
        DO J=1,NAT4
                    IF(IBASE4(J).EQ.INS(3)) NC8 = J
                    IF(IBASE4(J).EQ.INS(4)) NC6 = J
                    IF(IBASE4(J).EQ.IPS(3)) NN9 = J
                    IF(IBASE4(J).EQ.IPS(4)) NN1 = J
           KJ = IBASE4(J)
           ATOM1(J) = ATNM(KJ)
           DO KC = 1,3
              X1(KC,J) = CR(KJ,KC)
           END DO
C          N1(J) = 0
         IF(ATOM1(J).EQ.'C1'' ') THEN
                    DO KC=1,3
                       RC1P(KC) = X1(KC,J)
                    END DO
	 ELSE
           N1(J) = J
         END IF
        END DO
C
         DO KJ=1,3
           C8C1(KJ) = XC1P(KJ) - RC8(KJ)
           C6C1(KJ) = RC1P(KJ) - RC6(KJ)
         END DO
       SUMANG = 0.0
       SUMDS1 = 0.0
       SUMDS2 = 0.0
       DO KJ=1,3
         SUMDS1 = C8C1(KJ) * C8C1(KJ) + SUMDS1
         SUMDS2 = C6C1(KJ) * C6C1(KJ) + SUMDS2
       END DO
         SUMDS1 = SQRT(SUMDS1)
         SUMDS2 = SQRT(SUMDS2)
       DO KJ=1,3
         C8C1(KJ) = C8C1(KJ)/SUMDS1
         C6C1(KJ) = C6C1(KJ)/SUMDS2
	  BPDIR(KJ,3) = C8C1(KJ) * 1.0
	  BPDIR(KJ,4) = C6C1(KJ) * 1.0      
       END DO
       C8C6DS = 0.0
       DSN1C1 = 0.0
       DSN9C1 = 0.0
       XN1C1 = 0.0
       XN9C1 = 0.0
       C1C1DS = 0.0
       DO KJ=1,3
        C8C6DS = C8C6DS + (RC8(KJ) - RC6(KJ)) * (RC8(KJ) - RC6(KJ))
	C1C1DS = C1C1DS + (XC1P(KJ)-RC1P(KJ))*(XC1P(KJ)-RC1P(KJ))
	DSN9C1 = DSN9C1 + (X(KJ,NN9)-XC1P(KJ))*(X(KJ,NN9)-XC1P(KJ))
	DSN1C1 = DSN1C1 + (X1(KJ,NN1)-RC1P(KJ))*(X1(KJ,NN1)-RC1P(KJ))
	XN9C1 = XN9C1 + (X(KJ,NN9)-XC1P(KJ))*(RC1P(KJ)-XC1P(KJ))
	XN1C1 = XN1C1 + (X1(KJ,NN1)-RC1P(KJ))*(XC1P(KJ)-RC1P(KJ))
       END DO
       OPENDS(LNDEX+2) = SQRT(C8C6DS)
       OPENC1(LNDEX+2) = SQRT(C1C1DS)
       DSN1C1 = SQRT(DSN1C1)
       DSN9C1 = SQRT(DSN9C1)
       ANGBN9(LNDEX+2) = ACOS(XN9C1/(DSN9C1*OPENC1(LNDEX+2))) * CONV
       ANGBN1(LNDEX+2) = ACOS(XN1C1/(DSN1C1*OPENC1(LNDEX+2))) * CONV
C----------------------------------------------------------------------
C      CHECKING  WHETHER THE MOLECULE GOES UP OR DOWN IN Z.
C----------------------------------------------------------------------
        DO K=1,3
          ORIG2(K) = (RC8(K) + RC6(K)) * 0.5
          VEC(K) = ORIG2(K) - ORIG(K)
        END DO
C
      ELSE
C
       DOUBL2 = .FALSE.
       IF(NC8.EQ.0.AND.NC6.NE.0) THEN
           NC8 = NC6
           DO KL=1,3
              RC8(KL) = RC6(KL)
           END DO
       ENDIF
       DO KJ=1,3
         VEC(KJ) = RC8(KJ) - XC8(KJ)
       END DO
C
      ENDIF
C
        CALL NORMAL(VEC(1),VEC(2),VEC(3),EK1,EK2,EK3)
        ANG = EL * EK1 + EM * EK2 + EN * EK3
       IF(ANG.LT.0.0) THEN
             EL = -EL
             EM = -EM
             EN = -EN
          IF(DOUBL1) THEN
             DO KL=1,3
                    BNORM(KL,1) = -BNORM(KL,1)
                    BNORM(KL,2) = -BNORM(KL,2)
             END DO
          ELSE
                  DO KL=1,3
                     YAXIS(1,KL) = -YAXIS(1,KL)
                  END DO
         ENDIF
      ENDIF
C
      IF(.NOT.DOUBL1) THEN
C 
      THETA = 10.0601
      IF(NAT1.EQ.10) THETA = 12.6921
         ELS = EL * 1.0
         EMS = EM * 1.0
         ENS = EN * 1.0
         CALL ROTMAT(ELS,EMS,ENS,THETA,R)
         CALL ROTVEC(R,XAXIS(1,1),XAXIS(1,2),XAXIS(1,3),TEMPX)
         CALL ROTVEC(R,YAXIS(1,1),YAXIS(1,2),YAXIS(1,3),TEMPY)
         DO KJ=1,3
                    YAXIS(1,KJ) = TEMPY(KJ)
                    XAXIS(1,KJ) = TEMPX(KJ)
                    ORIGIN(1,KJ) = XC8(KJ) - YAXIS(1,KJ) * 4.9122
         END DO
      ELSE 
C
      CALL CROSS(YAXIS(1,1),YAXIS(1,2),YAXIS(1,3),EL,EM,EN,XAX(1),

     1 XAX(2),XAX(3))
      CALL NORMAL(XAX(1),XAX(2),XAX(3),XAXIS(1,1),XAXIS(1,2),XAXIS(1,3))
      CALL DIRECX(XAXIS,1,XC1P1,XC81,YAXIS,NDX)
       IF(NDX.EQ.1) THEN
           DO KJ=1,3
              BNORM(KJ,1) = -BNORM(KJ,1)
              BNORM(KJ,2) = -BNORM(KJ,2)
           END DO
C           TMPBAS = BASE(NOBASE)(1:1)
C           BASE(NOBASE)(1:1) = BASE(NOBASE)(3:3)
C           BASE(NOBASE)(3:3) = TMPBAS
        END IF
      END IF
C
C   END OF  FIRST BASE/base-pair CALC.
C-----------------------------------------------------------------------
      IF(DOUBL2) THEN
C
C   START OF DOUBLE-STRANDED SECOND base-pair CALC.
      DO I=1,3
        ORIG1(I) = 0.0
        ORIG2(I) = 0.0
	IF(ANSOR.EQ.'Y') THEN
          DO NN=1,NAT3
	    IF(ATOM(NN) .NE. 'C1'' ')ORIG1(I) = ORIG1(I) + X(I,NN)
 	  END DO
            ORIG1(I) = ORIG1(I)/(NAT3-1)
	  DO NN=1,NAT4
	   IF(ATOM1(NN) .NE.'C1'' ')ORIG2(I) = ORIG2(I) + X1(I,NN)
	  END DO
            ORIG2(I) = ORIG2(I)/(NAT4-1)
            ORIG(I) = (ORIG1(I) + ORIG2(I)) * 0.5
        ELSE
            ORIG(I) = (RC8(I) + RC6(I)) * 0.5
	END IF
          ORIGIN(2,I) = ORIG(I)
          PROJ(I) =  RC8(I) - RC6(I)
      END DO
C
	IF(ANSC1.EQ.'Y') THEN
	   DO KKK=1,3
	      PROJ(KKK) = XC1P(KKK) - RC1P(KKK)
	      C8C6(KKK) = RC8(KKK) - RC6(KKK)
	      XMID(KKK) = (XC1P(KKK) + RC1P(KKK)) * 0.5
	   END DO
	CALL NORMAL(C8C6(1),C8C6(2),C8C6(3),YC8C6(1),YC8C6(2),YC8C6(3))
	CALL NORMAL(PROJ(1),PROJ(2),PROJ(3),YAXIS(2,1),YAXIS(2,2),YAXIS(
     1  2,3))
	   CONSTN = (YAXIS(2,1)*(XMID(1)-RC8(1)) + YAXIS(2,2)*(XMID(2)-
     1       RC8(2)) + YAXIS(2,3)*(XMID(3)-RC8(3)))/(YAXIS(2,1)*YC8C6(1)
     2       +YAXIS(2,2)*YC8C6(2) + YAXIS(2,3)*YC8C6(3))
	   XYZMD1 = RC8(1) + YC8C6(1) * CONSTN
	   XYZMD2 = RC8(2) + YC8C6(2) * CONSTN
	   XYZMD3 = RC8(3) + YC8C6(3) * CONSTN
	ELSE
	CALL NORMAL(PROJ(1),PROJ(2),PROJ(3),YAXIS(2,1),YAXIS(2,2),YAXIS(
     1   2,3))
	END IF
C
	IF((ANSOR.EQ.'Y').OR.(ANSC1.NE.'Y')) THEN
	ORIGIN(2,1) = ORIG(1)
	ORIGIN(2,2) = ORIG(2)
	ORIGIN(2,3) = ORIG(3)
	ELSE
	   ORIGIN(2,1) = XYZMD1
	   ORIGIN(2,2) = XYZMD2
	   ORIGIN(2,3) = XYZMD3
	END IF
	WRITE(9,'('' BASE-PAIR ORIGIN-1'',3F10.4)')(ORIGIN(2,KJ),KJ=1,3)
C ----------------------------------------------------------------------
C     FINDING BASE & base-pair NORMALS OF THE TOP PLANE
C ----------------------------------------------------------------------
      CALL PLANE(X,NAT3,N,1,DVA,EL2,EM2,EN2,DETR,DVS3)
C
      BNORM(1,3) = EL2
      BNORM(2,3) = EM2
      BNORM(3,3) = EN2
      CALL PLANE(X1,NAT4,N1,1,DVA,EL1,EM1,EN1,DETR,DVS4)
      ANG = EL2 * EL1 + EM2 * EM1 + EN2 * EN1
      IF(ANG.LT.0.0) THEN
         EL1 = -EL1
         EM1 = -EM1
         EN1 = -EN1
      ENDIF
      BNORM(1,4) = EL1
      BNORM(2,4) = EM1
      BNORM(3,4) = EN1
C
      EL1 = EL1 + EL2
      EM1 = EM1 + EM2
      EN1 = EN1 + EN2
      CALL NORMAL(EL1,EM1,EN1,EL2,EM2,EN2)
C
      ELSE
C----------------------------------------------------------------------
      CALL SINGAXS(IPS,N,NAT3,NC8,X,2,XAXIS,YAXIS,EL2,EM2,EN2,IBASE3,S2)
C
      ENDIF
C
      ANG = EL * EL2 +  EM * EM2 +  EN * EN2
      IF(ANG.LT.0.0) THEN
               EL2 = -EL2
               EM2 = -EM2
               EN2 = -EN2
               IF(.NOT.DOUBL2) THEN
                    DO KL=1,3
                       YAXIS(2,KL) = -YAXIS(2,KL)
                    END DO
               ELSE
                 DO KL=1,3
                    BNORM(KL,3) = -BNORM(KL,3)
                    BNORM(KL,4) = -BNORM(KL,4)
                 END DO
               ENDIF
       ENDIF
C
      IF(.NOT.DOUBL2) THEN
C
        THETA = 10.0601
        IF(NAT3.EQ.10) THETA = 12.6921
         ELS = EL2 * 1.0
         EMS = EM2 * 1.0
         ENS = EN2 * 1.0
         CALL ROTMAT(ELS,EMS,ENS,THETA,R)
         CALL ROTVEC(R,XAXIS(2,1),XAXIS(2,2),XAXIS(2,3),TEMPX)
         CALL ROTVEC(R,YAXIS(2,1),YAXIS(2,2),YAXIS(2,3),TEMPY)
         DO KJ=1,3
                    XAXIS(2,KJ) = TEMPX(KJ)
                    YAXIS(2,KJ) = TEMPY(KJ)
                    ORIGIN(2,KJ) = RC8(KJ) - YAXIS(2,KJ) * 4.9122
         END DO
C
      ELSE
C
      CALL CROSS(YAXIS(2,1),YAXIS(2,2),YAXIS(2,3),EL2,EM2,EN2,
     1XAX(1),XAX(2),XAX(3))
      CALL NORMAL(XAX(1),XAX(2),XAX(3),XAXIS(2,1),XAXIS(2,2),XAXIS(2,3))
      CALL DIRECX(XAXIS,2,XC1P2,XC82,YAXIS,NDX)
        IF(NDX.EQ.1) THEN
C              TMPBAS = BASE(NOBASE+1)(1:1)
C              BASE(NOBASE+1)(1:1) = BASE(NOBASE+1)(3:3)
C              BASE(NOBASE+1)(3:3) = TMPBAS
           DO KJ=1,3
              BNORM(KJ,3) = -BNORM(KJ,3)
              BNORM(KJ,4) = -BNORM(KJ,4)
           END DO 
         END IF
       END IF
C
C     FINDING THE X-AXIS OF THE TOP BASE PAIR PLANE
C
       DO  I=1,3
        BM(I) = ORIGIN(2,I) - ORIGIN(1,I)
       END DO
C
        WRITE(9,*) 'FIRST base-pair NORMAL',EL,EM,EN
	ZAXIS(1,1) = EL
	ZAXIS(1,2) = EM
	ZAXIS(1,3) = EN
 	WRITE(9,*) 'SECOND base-pair NORMAL', EL2,EM2,EN2
 	ZAXIS(2,1) = EL2
	ZAXIS(2,2) = EM2
	ZAXIS(2,3) = EN2
        LNDEX=LNDEX + 1
        DO KK1=1,3
	  BPNORML(LNDEX,KK1) = ZAXIS(1,KK1)
	  BPNORML(LNDEX+1,KK1) = ZAXIS(2,KK1)
	END DO
C
      WRITE(9,*) '(ORIGIN) CENTER TO CENTER VECTOR _'
      WRITE(9,'(3F8.3)') (BM(I),I=1,3)
C     WRITE(9,'(3F8.2)')((ORIGIN(J,I),I=1,3),J=1,2)
C      
      IF(DOUBL1) THEN
C
         WRITE(9,*) ' (1st B.P.) BASE - NORMALS ---'
	 WRITE(9,'(3F10.6,'' DV ='',F10.7)') (BNORM(M,1),M=1,3),DVS1
	 WRITE(9,'(3F10.6,'' DV ='',F10.7)') (BNORM(M,2),M=1,3),DVS2
C
         CALL BUCK(YAXIS,1,BNORM,PR1,PR2)
         CALL BUCK(XAXIS,1,BNORM,BK1,BK2)
         PROTW(LNDEX) = PR1
         BUCAN(LNDEX) = -BK1
	 CALL BUCK(ZAXIS,1,BPDIR,PR1,PR2)
         OPENAN(LNDEX) = PR2 
	ELSE
	WRITE(9,'('' Base Normal'',3F10.6,'' DV ='',F10.7)') EL,EM,EN,S1
        END IF
C
       IF(DOUBL2) THEN
C
         WRITE(9,*) ' (2nd. B.P.) BASE - NORMALS---'
	WRITE(9,'(3F10.6,'' DV ='',F10.7)') (BNORM(M,3),M=1,3),DVS3
	WRITE(9,'(3F10.6,'' DV ='',F10.7)') (BNORM(M,4),M=1,3),DVS4
         CALL BUCK(YAXIS,2,BNORM,PR2,PR1)
         CALL BUCK(XAXIS,2,BNORM,BK2,BK1)
         PROTW(LNDEX+1) = PR2
         BUCAN(LNDEX+1) = -BK2
	 CALL BUCK(ZAXIS,2,BPDIR,PR1,PR2)
	 OPENAN(LNDEX+1) = PR2 
       ELSE
      WRITE(9,'('' Base Normal'',3F10.6,'' DV ='',F10.7)')EL2,EM2,EN2,S2
      END IF
C
      CALL WEGDM(XAXIS,YAXIS,BM,STKINV,NRTCAL)
      WRITE(9,*) 'XAXIS(1)',(XAXIS(1,I),I=1,3)
      WRITE(9,*) 'XAXIS(2)',(XAXIS(2,I),I=1,3)
      WRITE(9,*) 'YAXIS(1)',(YAXIS(1,I),I=1,3)
      WRITE(9,*) 'YAXIS(2)',(YAXIS(2,I),I=1,3)
      DO KJ=1,3
      BSCENTR(LNDEX,KJ) = ORIGIN(1,KJ)
      BSCENTR((LNDEX+1),KJ) = ORIGIN(2,KJ)
      HORIG(LNDEX,KJ) = BSCENTR(LNDEX,KJ) - DXLOC(LNDEX) * XAXIS(1,KJ)
     1                - DYLOC(LNDEX) * YAXIS(1,KJ)
      END DO
      IF(STKINV) CALL WEGDM(XAXIS,YAXIS,BM,STKINV,NRTCAL)
C      IF(NTIM .EQ. 1 .AND. (ANS .EQ. 'Y'))THEN
        CALL GLOBAL(ORIGIN,XAXIS,YAXIS,LNDEX)
C      END IF
C ----------------------------------------------------------------------
C      CLEARING THE ARRAYS CONTAINING COORDINATES TO 0.0 BEFORE GOING
C              TO THE NEXT base-pair STEP.
C ----------------------------------------------------------------------
      DO K=1,3
          XC8(K) = 0.0
          XC6(K) = 0.0
          SM(K) = 0.000000000000000
          RC8(K) = 0.0
          RC6(K) = 0.0
          BM(K) = 0.000000000000000
          PROJ(K) = 0.0
          ORIG(K) = 0.0
          BMN(K) = 0.000000000000000
          XYZ1(K) = 0.0
          XYZ2(K) = 0.0
C         AX(K) = 0.0 $ C(K) = 0.0
          XYZ3(K) = 0.0
          XYZ4(K) = 0.0
      END DO
      DO K=1,15
           NSP(K) = 0
           NSP1(K) = 0
           N(K) = 0
           N1(K) = 0
           DVA(K) = 0.0
                    DO J=1,3
                       X(J,K) = 0.0
                       X1(J,K) = 0.0
                    END DO
      END DO
	WRITE(9,'(3(''****************************''))')
      RETURN
      END
C ********************************************************************
 
      SUBROUTINE CROSS(A,B,C,A1,B1,C1,ALPHA,BETA,GAMA)
      DOUBLE PRECISION A,B,C,ALPHA,BETA,GAMA
      DOUBLE PRECISION A1,B1,C1
      ALPHA=B*C1-C*B1
      BETA=C*A1-A*C1
      GAMA=A*B1-B*A1
      AMOD=DSQRT(ALPHA*ALPHA+BETA*BETA+GAMA*GAMA)
      ALPHA=ALPHA/AMOD
      BETA=BETA/AMOD
      GAMA=GAMA/AMOD
      RETURN
      END
 
C --------------------------------------------------------------------C
C  NORMALISING SUBROUTINE
C --------------------------------------------------------------------
      SUBROUTINE NORMAL(A,B,C,X,Y,Z)
      DOUBLE PRECISION A,B,C
      DOUBLE PRECISION X,Y,Z
      DOUBLE PRECISION ALEN,ALENGT

      ALEN =  A * A + B * B + C * C

      ALENGT = DSQRT(ALEN)

	IF (ABS(ALENGT) .LT. 1.0E-6) THEN
	 X=0.D0
	 Y=0.D0
	 Z=0.D0
	ELSE
         X = A/ALENGT
         Y = B/ALENGT
         Z = C/ALENGT
	ENDIF

      RETURN
      END
C --------------------------------------------------------------------
C       SUBROUTINE FOR DETERMINATION OF BEST PLANE NORMAL FROM THE
C       NO. OF POINTS.(MODIFIED TO TAKE CARE OF PERFECT PLANES)!!
C --------------------------------------------------------------------
      SUBROUTINE PLANE(X,NTOM,NX,NPLN,DV,EL,EM,EN,DET,DVS)
      DOUBLE PRECISION EL,EM,EN
      DIMENSION X(3,99),NX(99),DV(99),Y(3,99),T(3,3),B(3),A(3)
      DIMENSION XY(3),XZ(3),VEC1(3),VEC2(3)
      INTEGER XY,XZ
        DMOVE = 0.0
C
	NCH = 0
 111 	CONTINUE
	NCOUNT = 0
C
       DO 100 J1=1,NTOM
         NX1=NX(J1)
	 IF(NX1.NE.0) THEN
            NCOUNT = NCOUNT + 1
	    DO  K=1,3
	       Y(K,NCOUNT)=X(K,NX1) + DMOVE
	    END DO
	 END IF
 100  CONTINUE
	NTOM1 = NCOUNT
      DO 101 J1=1,3
      B(J1)=0.0
      DO 101 J2=1,3
      T(J1,J2)=0.0
 101  CONTINUE
      DO  102 J1=1,NTOM1
      T(1,1)=T(1,1)+Y(1,J1)*Y(1,J1)
      T(1,2)=T(1,2)+Y(1,J1)*Y(2,J1)
      T(1,3)=T(1,3)+Y(1,J1)*Y(3,J1)
      T(2,2)=T(2,2)+Y(2,J1)*Y(2,J1)
      T(2,3)=T(2,3)+Y(2,J1)*Y(3,J1)
      T(3,3)=T(3,3)+Y(3,J1)*Y(3,J1)
      B(1)=B(1)+Y(1,J1)
      B(2)=B(2)+Y(2,J1)
      B(3)=B(3)+Y(3,J1)
 102  CONTINUE
      DO 103 J1=1,3
      DO 103 J2=J1,3
      T(J2,J1)=T(J1,J2)
 103  CONTINUE
C     WRITE(6,3) ((T(I,J),J=1,3),I=1,3)
C     WRITE(6,2) (B(I),I=1,3)
  3     FORMAT(3X,'COEFFICIENTS',3F14.4)
  2     FORMAT(/5X,'CONSTANTS',3F14.4)
      CALL MINV(T,3,DX,XY,XZ)
      IF(DX.NE.0.0E0) THEN
      B1 = B(1)
      B2 = B(2)
      B3 = B(3)
      CALL MATMUL(T,B1,B2,B3,A)
      DET=1.0/(A(1)*A(1)+A(2)*A(2)+A(3)*A(3))
      DET=SQRT(DET)
      EL=A(1)*DET
      EM=A(2)*DET
      EN=A(3)*DET
	DVS = 0.0
      DO 104 J1=1,NTOM1
      DV(J1)=EL*Y(1,J1)+EM*Y(2,J1)+EN*Y(3,J1)-DET
	DVS = DVS + DV(J1) * DV(J1)
 104  CONTINUE
	DVS = SQRT(DVS)
	IF(DVS.GT.0.75.AND.NPLN.NE.0) THEN
          DMOVE = DMOVE + 25.0
	GO TO 111
	END IF
C
      ELSE
      WRITE(6,*)' ALL ATOMS ARE IN-PLANE, MOLECULE SHIFTED FOR CALC.'
       DMOVE = 1.0 + DMOVE
      GO TO 111
C      STOP
C 
      ENDIF
      RETURN
      END
C --------------------------------------------------------------------
C    SUBROUTINE FOR TIP AND INCL. CALC. W.R.T. LOCAL HELIX AXIS & 
C     WTILT AND WROLL CALC. W.R.T. THE MEAN DOUBLET AXES. 
C --------------------------------------------------------------------
      SUBROUTINE WEGDM(XAXIS,YAXIS,BMN,STKINV,NRTCAL)
      DOUBLE PRECISION XAXIS,YAXIS,ZAXIS,XM,YM,ZM1,ZM2,ZM3,Y1,Y2,Y3,
     1   Y12,Y22,Y32,X11,X12,X13,BMN,TL2,RL2,WTILT,WROLL,Y11,Y13
      COMMON /LOCAL/TILTL(99),ROLLL(99),TWISTL(99),SLY(99),SLX(99),
     1DZL(99),LNDEX,DIETLT(99),DIEROL(99),PROTW(99),BUCAN(99),SLZ(99),
     2TWISTH(99),OPENAN(99),OPENDS(99),ANGBN9(99),OPENC1(99),ANGBN1(99)
      DOUBLE PRECISION ZAXISH
      COMMON /ZSTRAXS/ZAXISH(3,99)
      DIMENSION XAXIS(2,3),YAXIS(2,3),XM(3),YM(3),ZAXIS(3),BMN(3)
      COMMON /DISPL/DXLOC(99),DYLOC(99),DZLOC(99),KKTI,
     1       BSCENTR(99,3),HORIG(99,3)
	COMMON /BASINF/BASE(99)
      CHARACTER*3 BASE
      LOGICAL STKINV
      CONV=180.0/3.141592654
      XM(1) = XAXIS(1,1) - XAXIS(2,1)
      XM(2) = XAXIS(1,2) - XAXIS(2,2)
      XM(3) = XAXIS(1,3) - XAXIS(2,3)
      YM(1) = YAXIS(1,1) - YAXIS(2,1)
      YM(2) = YAXIS(1,2) - YAXIS(2,2)
      YM(3) = YAXIS(1,3) - YAXIS(2,3)
      CALL CROSS(XM(1),XM(2),XM(3),YM(1),YM(2),YM(3),ZM1,ZM2,ZM3)
      CALL NORMAL(ZM1,ZM2,ZM3,ZAXIS(1),ZAXIS(2),ZAXIS(3))
      DO KJ=1,3
      ZAXISH(KJ,LNDEX) = ZAXIS(KJ)
      END DO
      CALL CROSS(XAXIS(1,1),XAXIS(1,2),XAXIS(1,3),ZAXIS(1),ZAXIS(2),ZAXI
     1S(3),Y1,Y2,Y3)
      CALL CROSS(XAXIS(2,1),XAXIS(2,2),XAXIS(2,3),ZAXIS(1),ZAXIS(2),ZAXI
     1S(3),Y12,Y22,Y32)
      ANG = Y1 * Y12 + Y2 * Y22 + Y3 * Y32
      ANGLE = ACOS(ANG) * CONV
      IF(ABS(ANGLE).GT.90.0.AND.(.NOT.STKINV)) STKINV = .TRUE.
      CALL CROSS(Y1,Y2,Y3,Y12,Y22,Y32,ZM1,ZM2,ZM3)
      ANG = ZM1 * ZAXIS(1) + ZM2 * ZAXIS(2) + ZM3 * ZAXIS(3)
C   
      IF(ANG.LT.0) ANGLE = -ANGLE
      TWISTH(LNDEX) = ANGLE
C     WRITE(9,8) ANGLE
      WRITE(9,99)
 99   FORMAT(/)
      WRITE(9,*) 'Z*-AXIS',(ZAXIS(I),I=1,3)
       AINC = 0.0
       TIP1 = 0.0
         TSLIDE = 0.0
      DO I=1,3
            AINC = AINC + YAXIS(1,I) * ZAXIS(I)
            TIP1 = TIP1 + XAXIS(1,I) * ZAXIS(I)
            TSLIDE = TSLIDE + BMN(I) * ZAXIS(I)
      END DO
C  
      DZL(LNDEX) = TSLIDE
       AINCL =    ASIN(AINC) * CONV
        TIPL =  -(ASIN(TIP1) * CONV)
      DIETLT(LNDEX) = AINCL
        DIEROL(LNDEX) = TIPL
C 
      XM(1) = XAXIS(1,1) + XAXIS(2,1)
      XM(2) = XAXIS(1,2) + XAXIS(2,2)
      XM(3) = XAXIS(1,3) + XAXIS(2,3)
      YM(1) = YAXIS(1,1) + YAXIS(2,1)
      YM(2) = YAXIS(1,2) + YAXIS(2,2)
      YM(3) = YAXIS(1,3) + YAXIS(2,3)
      CALL CROSS(XM(1),XM(2),XM(3),YM(1),YM(2),YM(3),ZM1,ZM2,ZM3)
      CALL NORMAL(ZM1,ZM2,ZM3,ZAXIS(1),ZAXIS(2),ZAXIS(3))
      CALL NORMAL(XM(1),XM(2),XM(3),X11,X12,X13)
      CALL NORMAL(YM(1),YM(2),YM(3),Y11,Y12,Y13)
C
      XM(1) = X11
        XM(2) = X12
         XM(3) = X13
      YM(1) = Y11
        YM(2) = Y12
         YM(3) = Y13
      TL2 = 0.0
        RL2 = 0.0
         CNEW = 0.0
          TSLIDE = 0.0
           SLIDE = 0.0
      WRITE(9,*) 'MEAN Z-AXIS ',(ZAXIS(I),I=1,3)
      DO I=1,3
            TL2 = TL2 + YAXIS(1,I) * ZAXIS(I)
            RL2 = RL2 + XAXIS(1,I) * ZAXIS(I)
                    CNEW = CNEW + XM(I) * BMN(I)
                    TSLIDE = TSLIDE + BMN(I) * ZAXIS(I)
                    SLIDE = SLIDE + BMN(I) * YM(I) 
      END DO
C    
      WTILT = - 2.0 * DASIN(TL2) * CONV
      WROLL =  2.0 * DASIN(RL2) * CONV
C    
      TILTL(LNDEX) = WTILT
        ROLLL(LNDEX) = WROLL
      SLY(LNDEX) = SLIDE
        SLX(LNDEX) = CNEW
         SLZ(LNDEX) = TSLIDE
      CALL CROSS(XAXIS(1,1),XAXIS(1,2),XAXIS(1,3),ZAXIS(1),ZAXIS(2),ZAXI
     1S(3),Y1,Y2,Y3)
      CALL CROSS(XAXIS(2,1),XAXIS(2,2),XAXIS(2,3),ZAXIS(1),ZAXIS(2),ZAXI
     1S(3),Y12,Y22,Y32)
      ANG = Y1 * Y12 + Y2 * Y22 + Y3 * Y32
      ANGLE = ACOS(ANG) * CONV
C   
      CALL CROSS(Y1,Y2,Y3,Y12,Y22,Y32,ZM1,ZM2,ZM3)
      ANG = ZM1 * ZAXIS(1) + ZM2 * ZAXIS(2) + ZM3 * ZAXIS(3)
      IF(ANG.LT.0.0) ANGLE = - ANGLE
      TWISTL(LNDEX) = ANGLE
      CALL SLIDCAL
C
       IF(STKINV) THEN
        DO KKN = 1,3
          XAXIS(2,KKN) = -XAXIS(2,KKN)
          YAXIS(2,KKN) = -YAXIS(2,KKN)
        END DO
       END IF
	IF(NRTCAL.NE.0) STKINV = .FALSE.
	NRTCAL = NRTCAL + 1
      RETURN
      END
C ----------------------------------------------------------------------
C       CALCULATION OF PROPELLER TWIST AND BUCKLE FROM THE BASE-NORMALS.
C ----------------------------------------------------------------------
      SUBROUTINE BUCK(AXIS,I,BNORML,VALUE2,VALUE)
      DOUBLE PRECISION AXIS,PROP,BNORML,P1,P2,P4,P5
      DIMENSION AXIS(2,3),PROP(3,2),BNORML(3,4),P1(3),P2(3),P4(3),P5(3)
      CONV=180.0/3.141592654
      J=1
      IF(I.EQ.2) J=3
      J1 = J + 1
C
C  CALCULATION OF PARAMETERS TAKING THE PROJECTION BY METHOD OF
C  "GRAM-SMGITH".
C --------------------------------------------------------------------
       PROJEC = 0.0
       PROJE2 = 0.0
       PRTW = 0.0
       SIGN = 0.0
       DO KJ=1,3
       PROJEC = PROJEC + AXIS(I,KJ) * BNORML(KJ,J)
       PROJE2 = PROJE2 + AXIS(I,KJ) * BNORML(KJ,J1)
       END DO
       DO KJ=1,3
       PROP(KJ,1) = BNORML(KJ,J) - PROJEC * AXIS(I,KJ)
       PROP(KJ,2) = BNORML(KJ,J1) - PROJE2 * AXIS(I,KJ)
       END DO
       CALL NORMAL(PROP(1,1),PROP(2,1),PROP(3,1),P1(1),P1(2),P1(3))
       CALL NORMAL(PROP(1,2),PROP(2,2),PROP(3,2),P2(1),P2(2),P2(3)) 
       DO KJ=1,3
       PRTW = PRTW + P1(KJ) * P2(KJ)
       END DO
       VALUE = ACOSF(PRTW) 
C --------------------------------------------------------------------
C   CALCULATION OF PARAMETERS BY TAKING CROSS PRODUCT.
C---------------------------------------------------------------------
      CALL CROSS(AXIS(I,1),AXIS(I,2),AXIS(I,3),BNORML(1,J),BNORML(2,J),
     1BNORML(3,J),P4(1),P4(2),P4(3))
      CALL CROSS(AXIS(I,1),AXIS(I,2),AXIS(I,3),BNORML(1,J1),BNORML(2,J1)
     1,BNORML(3,J1),P5(1),P5(2),P5(3))
      PRTW = 0.0
      DO KJ=1,3
      PRTW = PRTW + P4(KJ) * P5(KJ)
      END DO
      VALUE2 = ACOSF(PRTW) 
      SIGN = 0.0
      CALL CROSS(P4(1),P4(2),P4(3),P5(1),P5(2),P5(3),P1(1),P1(2),P1(3))
      DO KJ=1,3
      SIGN = SIGN + P1(KJ) * AXIS(I,KJ)
      END DO
      IF(SIGN.GT.0.0) VALUE2 = - VALUE2
      RETURN
      END
C----------------------------------------------------------------------
C  SUBROUTINE TO FIND OUT DIRECTION TOWARD MAJOR GROOVE, TO FIX X-AXIS.
C----------------------------------------------------------------------
      SUBROUTINE DIRECX(AXIS,I,XC8,XN9,Y,IND)
      DOUBLE PRECISION AXIS,VEC,VECN,Y
      DIMENSION AXIS(2,3),XC8(3),XN9(3),VEC(3),VECN(3),Y(2,3)
        IND = 0
      DO JK = 1,3
      VEC(JK) = XC8(JK) - XN9(JK)
      END DO
      CALL NORMAL(VEC(1),VEC(2),VEC(3),VECN(1),VECN(2),VECN(3))
      ANG = 0.0
      DO JK = 1,3
      ANG = ANG + AXIS(I,JK) * VECN(JK)
      END DO
C
      IF(ANG.GT.0.0) THEN
C
        IND = 1
           DO KJ=1,3
                      AXIS(I,KJ) = -AXIS(I,KJ)
                      Y(I,KJ) = -Y(I,KJ)
           END DO
      ENDIF
      RETURN
      END
C ---------------------------------------------------------------------
C
C     SUBROUTINE FOR CALCULATION OF HELICAL (GLOBAL) PARAMETERS.
C
C ---------------------------------------------------------------------
        SUBROUTINE GLOBAL(ORIG,XAXIS,YAXIS,LNDEX)
        DOUBLE PRECISION XAXIS,YAXIS
        DIMENSION ORIG(2,3),XAXIS(2,3),YAXIS(2,3)
      COMMON /GLOB/TILTG(99),ROLLG(99),TWISTG(99),DXG(99),DYG(99),HG(99)
        CONV = 180.0/3.1415926
C
        HREPIT = ORIG(1,3) - ORIG(2,3)
        IF(HREPIT.LT.0.0) HREPIT = -HREPIT
        HG(LNDEX) = HREPIT
C
        CALL DXDY(ORIG,YAXIS,1,DX1,DY1)
        CALL DXDY(ORIG,YAXIS,2,DX2,DY2)
        DXG(LNDEX) = DX1
        DYG(LNDEX) = DY1
        DXG(LNDEX+1) = DX2
        DYG(LNDEX+1) = DY2
C
         TILT = ASIN(YAXIS(1,3)) * CONV
         ROLL = ASIN(-XAXIS(1,3)) * CONV
         TILT2 = ASIN(YAXIS(2,3)) * CONV
         ROLL2 = ASIN(-XAXIS(2,3)) * CONV
       TILTG(LNDEX) = TILT
        ROLLG(LNDEX) = ROLL
       TILTG(LNDEX+1) = TILT2
        ROLLG(LNDEX+1) = ROLL2
C
         DIM = YAXIS(1,1) *YAXIS(1,1) + YAXIS(1,2) * YAXIS(1,2)
         DIM = SQRT(DIM)
          Y1Z1 = YAXIS(1,1) / DIM
          Y1Z2 = YAXIS(1,2) / DIM
C
         DIM = YAXIS(2,1) * YAXIS(2,1) + YAXIS(2,2) * YAXIS(2,2)
         DIM = SQRT(DIM)
          Y2Z1 = YAXIS(2,1) / DIM
          Y2Z2 = YAXIS(2,2) / DIM
          ANG = Y1Z1 * Y2Z1 + Y1Z2 * Y2Z2
          ANGLE = ACOS(ANG) * CONV
         TWISTG(LNDEX) = ANGLE
	 ZCOM = Y1Z1 * Y2Z2 - Y1Z2 * Y2Z1
	 IF(ZCOM.LT.0.0) TWISTG(LNDEX) = -TWISTG(LNDEX)
        RETURN
        END
C
        SUBROUTINE DXDY(ORIGIN,AXIS,I,DX,DY)
        DOUBLE PRECISION AXIS
C       REAL LY2MY2,NUM
        DIMENSION ORIGIN(2,3),AXIS(2,3)
C
         Y2LY2M = AXIS(I,1) * AXIS(I,1) + AXIS(I,2) * AXIS(I,2)
        GRAD = AXIS(I,2) / AXIS(I,1)
        ANUM = ORIGIN(I,2) - GRAD * ORIGIN(I,1)
C
        DX = -AXIS(I,1) * ANUM /(SQRT(Y2LY2M))
C 
          X1 = -ANUM * AXIS(I,1) * AXIS(I,2) / Y2LY2M
          Y1 = AXIS(I,1) * AXIS(I,1) * ANUM / Y2LY2M
C    
          DY = (X1-ORIGIN(I,1)) * (X1-ORIGIN(I,1)) + (Y1-ORIGIN(I,2)) *
     1 (Y1-ORIGIN(I,2))
          DY = SQRT(DY)
        EL = (X1 - ORIGIN(I,1))/DY
        EM = (Y1 - ORIGIN(I,2))/DY
      ANG = EL * AXIS(I,1) + EM * AXIS(I,2)
      IF(ANG.GT.0.0) DY = -DY
        RETURN
        END
C
C-----------------------------------------------------------------------
      SUBROUTINE WRITING(FILENM,TITLE1,INPUT,NTIM,NRUN,ANORIN,ANDPDB)
      COMMON /GLOB/TILTG(99),ROLLG(99),TWISTG(99),DXG(99),DYG(99),HG(99)
      COMMON /LOCAL/TILTL(99),ROLLL(99),TWISTL(99),SLY(99),SLX(99),
     1DZL(99),LNDEX,DIETLT(99),DIEROL(99),PROTW(99),BUCAN(99),SLZ(99),
     2TWISTH(99),OPENAN(99),OPENDS(99),ANGBN9(99),OPENC1(99),ANGBN1(99)
      DOUBLE PRECISION XAXIS,YAXIS,ZAXIS,XM,YM,ZM1,ZM2,ZM3,Y1,Y2,Y3,
     1   Y12,Y22,Y32,X11,X12,X13,BMN,TL2,RL2,WTILT,WROLL,Y11,Y13
      COMMON /DISPL/DXLOC(99),DYLOC(99),DZLOC(99),KKTI,
     1       BSCENTR(99,3),HORIG(99,3)
	COMMON /BASINF/BASE(99)
      DOUBLE PRECISION ZAXISH,BPNORML
      COMMON /BPZAXS/BPNORML(99,3)
      COMMON /ZSTRAXS/ZAXISH(3,99)
      COMMON /GLOREO/DCRS1,DCRS2,DCRS3,ANGLE,DX1,DY1,INFBP(2,99)
      COMMON /REPLY/ANSOR,ANSC1,ANSBPN,ANSDB,ANS,ANSHO
      DIMENSION ANGST(99,99),TITLE1(20)
      CHARACTER *80 TITLE,FILENM,INPUT*40
      CHARACTER *4 ANORIN
      CHARACTER*3 BASE
      CHARACTER*1 ANSOR,ANSC1,ANS,ANSHO,ANSDB,ANSBPN
      CHARACTER*80 TITLE1
      CONV=180.0/3.141592654

      IF(NTIM .EQ. 1)THEN
C
	KKKK=INDEX(FILENM(1:80),']')

	IF (KKKK .EQ. 0) THEN
	 KKKK=1
	ELSE
	 KKKK=KKKK+1
	ENDIF

	KKKKK=INDEX(FILENM(KKKK:80),'.')
	KKKKK=KKKKK+KKKK

	TITLE=FILENM(KKKK:KKKKK-1)//'PRM'

      OPEN(UNIT=11,FILE=TITLE,STATUS='NEW')
      WRITE(11,300) FILENM
 300  FORMAT(1X,'COORDINATES INPUT TAKEN FROM ::',A80)
 400  FORMAT(/' Results of the calculation are written into: '/1X,A80)
          WRITE(6,400) TITLE
	  WRITE(8,400) TITLE
       DO II=1,KKTI
          WRITE(11,'(A80)') TITLE1(II)
       END DO
        IF(ANSBPN.NE.'Y') WRITE(11,61)INPUT
 61   FORMAT(/1X,'RESIDUE NUMBERS TAKEN FROM ::',A40)

	IF(ANSOR.EQ.'Y') WRITE(11,71)
 71   FORMAT(/2X,'*Basepair origin is at the center of gravity*')

	IF(ANSC1.EQ.'Y') THEN
          WRITE(11,81)
 81   FORMAT(/2X,'The basepair Y-axis is along C1''--C1'' direction.'/)
        ELSE
          WRITE(11,82)
 82   FORMAT(/2X,'The basepair Y-axis is along C6--C8 direction.'/)
        END IF
      END IF
C-----------------------------------------------------------------------
      IF (NTIM .EQ. 2)WRITE(11,221)
 221  FORMAT('1',2X,'Single strand parameters for Strand 1 of duplex:')
      IF (NTIM .EQ. 3)WRITE(11,222)
 222  FORMAT('1',2X,'Single strand parameters for Strand 2 of duplex:')
      WRITE(11,1)
      WRITE(11,199)
      DO K=1,LNDEX
      WRITE(11,100) INFBP(1,K),BASE(K),INFBP(2,K), TILTL(K),ROLLL(K),
     1        TWISTL(K),SLX(K),SLY(K) ,SLZ(K)
      END DO
      CALL STATIS(TILTL,LNDEX,AVTL,STDTL)
      CALL STATIS(ROLLL,LNDEX,AVRL,STDRL)
      CALL STATIS(TWISTL,LNDEX,AVTW,STDTW)
      CALL STATIS(SLX,LNDEX,AVSX,STDSX)
      CALL STATIS(SLY,LNDEX,AVSY,STDSY)
      CALL STATIS(SLZ,LNDEX,AVSZ,STDSZ)
      CALL STATIS(PROTW,(LNDEX+1),AVPR,STDPR)
      CALL STATIS(BUCAN,(LNDEX+1),AVBU,STDBU)
      CALL STATIS(OPENAN,(LNDEX+1),AVAN,STDAN)
      CALL STATIS(OPENDS,(LNDEX+1),AVDS,STDDS)
      CALL STATIS(ANGBN9,(LNDEX+1),AVBN9,STDBN9) 
      CALL STATIS(ANGBN1,(LNDEX+1),AVBN1,STDBN1) 
      CALL STATIS(OPENC1,(LNDEX+1),AVBDC1,STDC1)
      WRITE(11,100) INFBP(1,(LNDEX+1)),BASE(LNDEX+1),INFBP(2,(LNDEX+1))
       WRITE(11,199)
       WRITE(11,110) AVTL,AVRL,AVTW,AVSX,AVSY,AVSZ
       WRITE(11,111) STDTL,STDRL,STDTW,STDSX,STDSY,STDSZ
       WRITE(11,199)
 199  FORMAT(' ---------------------------------------------------------
     1------------')
 999   FORMAT('1')
 110   FORMAT(5X,'Av.',4X,10F8.2)
 111   FORMAT(3X,'Std.Dev. ',10F8.2)
  1   FORMAT(2X,'Local Step Parameters: '/16X,'Tilt    Roll   Twist     
     1Dx      Dy      Dz')
  2   FORMAT(/2X,'GLOBAL Helical Parameters:',/14X,'Inclin.   Tip    Twi    
     1st','     dx      dy      dz')
  3   FORMAT(/,2X,'Local Helical Parameters:',/14X,'Inclin.   Tip    Twi
     1st',5X,'dx',6X,'dy',6X,'dz')
 100  FORMAT(1X,I2,2X,A3,2X,I2,10(1X,F7.2))
 101  FORMAT(1X,I2,2X,A3,2X,I2,2F8.2,8X,2F8.2)
      WRITE(11,3)
      WRITE(11,199)
      DO K=1,LNDEX
      WRITE(11,100) INFBP(1,K),BASE(K),INFBP(2,K),DIETLT(K),DIEROL(K),
     1             TWISTH(K),DXLOC(K),DYLOC(K),DZLOC(K)
      END DO
      WRITE(11,100) INFBP(1,LNDEX+1),BASE(LNDEX+1),INFBP(2,LNDEX+1)
      CALL STATIS(DIETLT,LNDEX,AVTL,STDTL)
      CALL STATIS(DIEROL,LNDEX,AVRL,STDRL)
      CALL STATIS(TWISTH,LNDEX,AVTW,STDTW)
      CALL STATIS(DXLOC,LNDEX,AVSX,STDSX)
      CALL STATIS(DYLOC,LNDEX,AVSY,STDSY)
      CALL STATIS(DZLOC,LNDEX,AVSZ,STDSZ)
      WRITE(11,199)
      WRITE(11,110) AVTL,AVRL,AVTW,AVSX,AVSY,AVSZ
      WRITE(11,111) STDTL,STDRL,STDTW,STDSX,STDSY,STDSZ
      WRITE(11,199)
C-----------------------------------------------------------------------
	IF(NTIM.LE.1) THEN
          IF(ANS .NE.'Y')THEN
          WRITE(11,90)
  	ELSE
 90   FORMAT('1  Molecule has NOT been reoriented along GLOBAL axis.')
       IF(ANORIN .EQ. 'C   ') THEN 
          WRITE(11,92)
 92	FORMAT('1  Molecule has been reoriented along best linear axis
     1 obtained from '/'  basepair centers.')
       ELSE IF(ANORIN.EQ.'O   ') THEN
          WRITE(11,91)
 91     FORMAT('1  Molecule has been reoriented along best linear axis
     1 obtained from '/'  local helix ORIGINs.')
       ELSE
	   WRITE(11,192)ANORIN
 192	FORMAT('1  Molecule has been reoriented along best linear axis
     1 obtained from ',A4,' atoms.')
        END IF
	WRITE(11,93) DCRS1,DCRS2,DCRS3,ANGLE
 	WRITE(11,94) DX1,DY1
       END IF
C
      END IF
900    CONTINUE
        WRITE(11,2)
        WRITE(11,199)
 93   FORMAT(/2X,'Molecule rotated along(',3F8.4,') by',F8.3,' Deg.')
 94   FORMAT(2X,'and translated by',2F7.3,'A along X and Y directions.')
      DO K=1,LNDEX
      WRITE(11,100)INFBP(1,K), BASE(K),INFBP(2,K),TILTG(K),ROLLG(K),
     1             TWISTG(K),DXG(K),DYG(K),HG(K)
      END DO
      KK = LNDEX + 1
      CALL STATIS(TILTG,(LNDEX+1),AVTL,STDTL)
      CALL STATIS(ROLLG,(LNDEX+1),AVRL,STDRL)
      CALL STATIS(TWISTG,(LNDEX),AVTW,STDTW)
      CALL STATIS(DXG,(LNDEX+1),AVSX,STDSX)
      CALL STATIS(DYG,(LNDEX+1),AVSY,STDSY)
      CALL STATIS(HG,(LNDEX),AVSZ,STDSZ)
      WRITE(11,101) INFBP(1,KK),BASE(KK),INFBP(2,KK),TILTG(KK),ROLLG(KK)
     1,DXG(KK),DYG(KK)
      WRITE(11,199)
      WRITE(11,110) AVTL,AVRL,AVTW,AVSX,AVSY,AVSZ
      WRITE(11,111) STDTL,STDRL,STDTW,STDSX,STDSY,STDSZ
      WRITE(11,199)
 950  CONTINUE
C-----------------------------------------------------------------------
  5     FORMAT(/2X,'Intra base-pair Parameters:'//15X,'PROP.  BUCKLE  OP
     1ENING   GLYCOSIDIC     C8..C6  C1..C1')
C-- NOT PRINTED OUT FOR SINGLE STRAND CALCULATION.
       IF((ANSDB.EQ.'Y').AND.NTIM.LE.1) THEN
        WRITE(11,5)
        WRITE(11,199)
         DO K=1,(LNDEX+1)
       WRITE(11,100)INFBP(1,K),BASE(K),INFBP(2,K),PROTW(K),BUCAN(K),
     1            OPENAN(K),ANGBN9(K),ANGBN1(K),OPENDS(K),OPENC1(K)
         END DO
         WRITE(11,199)
         WRITE(11,110) AVPR,AVBU,AVAN,AVBN9,AVBN1,AVDS,AVBDC1
         WRITE(11,111) STDPR,STDBU,STDAN,STDBN9,STDBN1,STDDS,STDC1
         WRITE(11,199)
       END IF
 1000 CONTINUE
C-----------------------------------------------------------------------
      IF (NTIM .EQ. 1)WRITE(11,4)
  4   FORMAT('1',6X,' LOCAL HELIX ORIGINS',15X,'LOCAL HELIX AXES'/)
       DO K=1,LNDEX
         ZS1 = ZAXISH(1,K) * 1.0
         ZS2 = ZAXISH(2,K) * 1.0
         ZS3 = ZAXISH(3,K) * 1.0
       IF(NTIM .EQ. 1)WRITE(11,200)(HORIG(K,M),M=1,3),ZS1,ZS2,ZS3
       END DO
  201  FORMAT(51X,3F8.3)
C
       IF(NTIM .EQ. 1)WRITE(11,96)
  96   FORMAT(//6X,' BASE/B.P.  CENTERS',16X,'BASE/B.P. NORMALS'/)
	DO K=1,LNDEX+1
          ZS1 = BPNORML(K,1) * 1.0
          ZS2 = BPNORML(K,2) * 1.0
          ZS3 = BPNORML(K,3) * 1.0
	IF(NTIM .EQ.1)WRITE(11,200)(BSCENTR(K,M),M=1,3),ZS1,ZS2,ZS3
	END DO
 200  FORMAT(3X,3F8.3,5X,3F10.4)
      INDM = LNDEX - 1
      DO K=1,LNDEX 
	 DO KK=1,LNDEX
           SUM = 0.0
           DO KJ=1,3
              SUM = SUM + ZAXISH(KJ,K) * ZAXISH(KJ,KK)
           END DO
           ANG = ACOS(SUM) * CONV
           ANGST(K,KK) = ANG
C         WRITE(11,108) K,(K+1),ANG
         END DO
       END DO
        WRITE(11,999)
 	WRITE(11,97)
	DO KKK=1,LNDEX
	   WRITE(11,11) KKK,(ANGST(KKK,K),K=1,LNDEX)
	END DO
  97	FORMAT(3X,'ANGLE BETWEEN SUCCESSIVE LOCAL HELIX AXES:'/)
  11	FORMAT(2X,I2,'.',10F7.2,9(/5X,10F7.2))
C
         DO K=1,LNDEX + 1
	   DO KK=1,LNDEX + 1
	     SUM = 0.0
	     DO KJ = 1,3
	       SUM = SUM + BPNORML(K,KJ)*BPNORML(KK,KJ)
	     END DO
	     ANGST(K,KK) = ACOSF(SUM)
	   END DO
	END DO
  98	FORMAT(/3X,'ANGLE BETWEEN SUCCESSIVE BASE/BASE-PAIR NORMALS:'/)
	WRITE(11,98)
	DO KKK=1,LNDEX + 1
	   WRITE(11,11) KKK,(ANGST(KKK,K),K=1,LNDEX+1)
        END DO
      RETURN
      END
C----------------------------------------------------------------------
      SUBROUTINE SLIDCAL
      COMMON /LOCAL/TILTL(99),ROLLL(99),TWISTL(99),SLY(99),SLX(99),
     1DZL(99),LNDEX,DIETLT(99),DIEROL(99),PROTW(99),BUCAN(99),SLZ(99),
     2TWISTH(99),OPENAN(99),OPENDS(99),ANGBN9(99),OPENC1(99),ANGBN1(99)
      DOUBLE PRECISION XAXIS,YAXIS,ZAXIS,XM,YM,ZM1,ZM2,ZM3,Y1,Y2,Y3,
     1   Y12,Y22,Y32,X11,X12,X13,BMN,TL2,RL2,WTILT,WROLL,Y11,Y13
      COMMON /DISPL/DXLOC(99),DYLOC(99),DZLOC(99),KKTI,
     1       BSCENTR(99,3),HORIG(99,3)
	COMMON /BASINF/BASE(99)
      DIMENSION N(3),M(3),A(3,3)
      CHARACTER*3 BASE

      SLIDEX = SLX(LNDEX)
       SLIDEY = SLY(LNDEX)
        SLIDEZ = SLZ(LNDEX)
      TILT = DIETLT(LNDEX)
       ROLL = -DIEROL(LNDEX)
        TWIST = TWISTL(LNDEX)
      CONV = 3.141592654/180.000
      ROLL = ROLL * CONV
        TILT = TILT * CONV
         TWIST = TWIST * CONV
      CY = COS(ROLL)
        CX = COS(TILT)
      SY = SIN(ROLL)
        SX = SIN(TILT)
      CT = COS(TWIST)
        ST = SIN(TWIST)
      A(1,1) = 2.0 * CX * ST
      A(1,2) = 0.0
      A(1,3) = 2.0 * SX
      A(2,1) = 0.0
      A(2,2) = -2.0 * CY * ST
      A(2,3) = 2.0 * SY
      A(3,1) = -2.0 * CY * SX * ST
      A(3,2) = 2.0 * CX * SY * ST
      A(3,3) = CX * CY * (1.0 + CT)
      CALL MINV(A,3,DXX,N,M)
      DX = 1.0 + CT + SX * SX * ( 1.0 - CT)
      DX = SQRT(2.0 * DX)
      DY = 1.0 + CT + SY * SY * (1.0 - CT)
      DY = SQRT(2.0 * DY)
      B1 = SLIDEY * DY
      B2 = SLIDEX * DX
      B3 = (SLIDEZ * DX * DY)/2
      D1 = A(1,1) * B1 + A(1,2) * B2 + A(1,3) * B3
      D2 = A(2,1) * B1 + A(2,2) * B2 + A(2,3) * B3
      D3 = A(3,1) * B1 + A(3,2) * B2 + A(3,3) * B3
      DXLOC(LNDEX) = D1
        DYLOC(LNDEX) = D2
         DZLOC(LNDEX) = D3
      RETURN
      END
C
C***********************************************************************
      SUBROUTINE MINV (A,N,D,L,M)
C***********************************************************************
C     ----- STANDARD IBM MATRIX INVERSION ROUTINE -----
C
      DIMENSION A(1),L(1),M(1)
C
C     ----- SEARCH FOR LARGEST ELEMENT -----
C
      D = 1.E0
      NK = -N
      DO 80 K = 1,N
      NK = NK+N
      L(K) = K
      M(K) = K
      KK = NK+K
      BIGA = A(KK)
      DO 20 J = K,N
      IZ = N*(J-1)
      DO 20 I = K,N
      IJ = IZ+I
   10 IF( ABS(BIGA)- ABS(A(IJ))) 15,20,20
   15 BIGA = A(IJ)
      L(K) = I
      M(K) = J
   20 CONTINUE
C
C     ----- INTERCHANGE ROWS -----
C
      J = L(K)
      IF(J-K) 35,35,25
   25 KI = K-N
      DO 30 I = 1,N
      KI = KI+N
      HOLD = -A(KI)
      JI = KI-K+J
      A(KI) = A(JI)
   30 A(JI) = HOLD
C
C     ----- INTERCHANGE COLUMNS -----
C
   35 I = M(K)
      IF(I-K) 45,45,38
   38 JP = N*(I-1)
      DO 40 J = 1,N
      JK = NK+J
      JI = JP+J
      HOLD = -A(JK)
      A(JK) = A(JI)
   40 A(JI) = HOLD
C
C     ----- DIVIDE COLUMN BY MINUS PIVOT -----
C
   45 IF(BIGA) 48,46,48
   46 D = 0.E0
      GO TO 150
   48 DO 55 I = 1,N
      IF(I-K) 50,55,50

   50 IK = NK+I
      A(IK) = A(IK)/(-BIGA)
   55 CONTINUE
C
C     ----- REDUCE MATRIX -----
C
      DO 65 I = 1,N
      IK = NK+I
      HOLD = A(IK)
      IJ = I-N
      DO 65 J = 1,N
      IJ = IJ+N
      IF(I-K) 60,65,60
   60 IF(J-K) 62,65,62
   62 KJ = IJ-I+K
      A(IJ) = HOLD*A(KJ)+A(IJ)
   65 CONTINUE
C
C     ----- DIVIDE ROW BY PIVOT -----
C
      KJ = K-N
      DO 75 J = 1,N
      KJ = KJ+N
      IF(J-K) 70,75,70
   70 A(KJ) = A(KJ)/BIGA
   75 CONTINUE
C
C     ----- PRODUCT OF PIVOTS -----
C
      D = D*BIGA
C
C     ----- REPLACE PIVOT BY RECIPROCAL -----
C
      A(KK) = 1.E0/BIGA
   80 CONTINUE
C
C     ----- FINAL ROW AND COLUMN INTERCHANGE -----
C
      K = N
  100 K = (K-1)
      IF(K) 150,150,105
  105 I = L(K)
      IF(I-K) 120,120,108
  108 JQ = N*(K-1)
      JR = N*(I-1)
      DO 110 J = 1,N
      JK = JQ+J
      HOLD = A(JK)
      JI = JR+J
      A(JK) = -A(JI)
  110 A(JI) = HOLD
  120 J = M(K)
      IF(J-K) 100,100,125
  125 KI = K-N
      DO 130 I = 1,N
      KI = KI+N
      HOLD = A(KI)
      JI = KI-K+J
      A(KI) = -A(JI)
  130 A(JI) = HOLD
      GOTO 100
  150 RETURN
      END
C
      SUBROUTINE SINGAXS(IP,N,NATL,NC8,X,L,XAXIS,Y,EL,EM,EN,IBASE,DVSX)
         COMMON/ATINF/CR(2400,3),SECSQ(300),
     1             ISTS(300),IENS(300),NSEG,NAT
	COMMON /ATINFC/ATNM(2400),SECNM(300)
      DOUBLE PRECISION XAX,XAXIS,EL,EM,EN,Y,TEMPX,TEMPY
      DIMENSION IP(4),IN(4),XAX(3),YAXIS(2,3),ORIG(3),X(3,15),Y(2,3),
     1          XAXIS(2,3),N(30),DVA(30),TEMPX(3),TEMPY(3),IBASE(10)
     2         ,R(3,3)
      CHARACTER*3 SECNM
	CHARACTER*4 ATNM
      L2 = L * 2
      NADD = 0
      DO I=1,NATL
          IF(IBASE(I).EQ.IP(L2)) NC1P = I
      END DO
      DO KJ=1,3
         XAX(KJ) = X(KJ,NC8) - X(KJ,NC1P)
      END DO
      CALL NORMAL(XAX(1),XAX(2),XAX(3),XAXIS(L,1),XAXIS(L,2),XAXIS(L,3))
      CALL PLANE(X,NATL,N,1,DVA,EL,EM,EN,DETR,DVSX)
      CALL CROSS(EL,EM,EN,XAXIS(L,1),XAXIS(L,2),XAXIS(L,3),Y(L,1),Y(L,2)
     1          ,Y(L,3))
      RETURN
      END
C
      SUBROUTINE ROTVEC(RMAT,A1,A2,A3,B)
      DOUBLE PRECISION A1,A2,A3,B,DRMAT(3,3)
      DIMENSION RMAT(3,3),B(3)

C	CONVERT MATRIX TO DOUBLE PRECISION

	DO I=1,3
	  DO J=1,3
	   DRMAT(I,J)=DBLE(RMAT(I,J))
	  ENDDO
	ENDDO	
	
      DO I=1,3
      B(I) = DRMAT(I,1) * A1 + DRMAT(I,2) * A2 + DRMAT(I,3) * A3
      END DO

      RETURN
      END
C
      SUBROUTINE STATIS(VAR,NVAR,AVG,STD)
      DIMENSION VAR(1)
	IF(NVAR.GE.2) THEN
      SUM = 0.0
	NPOINT = 0
      DO 10 IV = 1,NVAR
      IF(VAR(IV).NE.0.000000) THEN
      SUM = SUM + VAR(IV)
      NPOINT = NPOINT + 1
      END IF
10    CONTINUE
      AVG = SUM/NPOINT
      SUM = 0.0
      DO 20 IV = 1, NVAR
      IF(VAR(IV).NE.0.0000000) THEN
      SUM = SUM + VAR(IV)**2
      END IF
20    CONTINUE
      SUM = SUM/NPOINT
      STD = SUM-AVG**2
      IF(STD .LT. 0.0) STD = 0.0
      IF(NPOINT.NE.0) STD = SQRTF(STD)
	ELSE
	   AVG = VAR(1)
	   STD = 0.0
	END IF
      RETURN
      END
C
      SUBROUTINE ROTMAT(EL,EM,EN,THETA,R)
      DIMENSION R(3,3)
      CON=3.14195/180.0
      THETR=THETA*CON
      CST= COS(THETR)
      SNT= SIN(THETR)
      CSTF=1.0-CST
      R(1,1)=EL*EL*CSTF+CST
      R(2,2)=EM*EM*CSTF+CST
      R(3,3)=EN*EN*CSTF+CST
      AA=EL*EM*CSTF
      BB=EN*SNT
      R(1,2)=AA-BB
      R(2,1)=AA+BB
      AA=EL*EN*CSTF
      BB=EM*SNT
      R(1,3)=AA+BB
      R(3,1)=AA-BB
      AA=EM*EN*CSTF
      BB=EL*SNT
      R(2,3)=AA-BB
      R(3,2)=AA+BB
      RETURN
      END
C
        SUBROUTINE MATMUL(R,A1,A2,A3,B)
        DIMENSION B(3),R(3,3)
        DO I=1,3
            B(I) = R(I,1) * A1 + R(I,2) * A2 + R(I,3) * A3
        END DO
        RETURN
        END
C
       SUBROUTINE POINT(NPTS,DX,DY)
	COMMON /CENTER/ CENTRS(99,3)
	DX = 0.0
	DY = 0.0
	DO KJ=1,NPTS
	   DX = DX + CENTRS(KJ,1)
	   DY = DY + CENTRS(KJ,2)
	END DO
	DX = DX/NPTS
	DY = DY/NPTS
	RETURN
	END
C
      SUBROUTINE REWRIT(FILENM,TITLE,KKTI)
         COMMON/ATINF/CR(2400,3),SECSQ(300),
     1             ISTS(300),IENS(300),NSEG,NAT
	COMMON /ATINFC/ATNM(2400),SECNM(300)
      CHARACTER*80 TITLE(20)
      CHARACTER*3 SECNM
      CHARACTER*4 ATNM
      CHARACTER*80 FILENM,FILNEW

	KKKK=INDEX(FILENM(1:80),']')

	IF (KKKK .EQ. 0) THEN
	 KKKK=1
	ELSE
	 KKKK=KKKK+1
	ENDIF

	KKKKK=INDEX(FILENM(KKKK:80),'.')
  	KKKKK=KKKKK+KKKK

	FILNEW=FILENM(KKKK:KKKKK-1)//'COOR'
C
      OPEN(UNIT=7,FILE=FILNEW,STATUS='NEW')
C 
      DO I=1,KKTI
               WRITE(7,'(A80)') TITLE(I)
      END DO
      K1 = 1
      NEND = IENS(K1)
      DO K=1,NAT
             IF(K.GT.NEND) THEN
                   K1 = K1 + 1
                   NEND = IENS(K1)
             END IF
             WRITE(7,101) ATNM(K),SECNM(K1),K1,(CR(K,M),M=1,3)
      END DO
  101 FORMAT('ATOM',8X,A4,1X,A3,2X,I4,4X,3F8.3)
       WRITE(6,10)FILNEW
       WRITE(8,10)FILNEW
  10  FORMAT(1X,'Reoriented atomic coordinates written into:'/1X,A80)
       CLOSE(7)
       RETURN
       END
C-----------------------------------------------------------------------
	SUBROUTINE TORSION(ANSDB,ANSPDB)
         COMMON/ATINF/CR(2400,3),SECSQ(300),
     1             ISTS(300),IENS(300),NSEG,NAT
C
	COMMON /ATINFC/ATNM(2400),SECNM(300)
        DIMENSION COOR1(3),COOR2(3),N(5),COOR3(3),NN(2400,4),ALPHA(99),
     1    S(200),XYZ1(3),XYZ2(3),XYZ3(3),XYZ4(3),BETA(99),GAMA(99)
     2    ,DELTA(99),EPSILN(99),ZETA(99),CHI(99),ANEU0(99),
     3   ANEU1(99),ANEU2(99),ANEU3(99),ANEU4(99),ANEUMAX(99),PHASE(99)
        DIMENSION XYZ(2400,3),EE(21),ISTSN(99),IENSN(99),SECNMN(99)
        CHARACTER*12 EE,CLASS
	CHARACTER *1 ANSDB,ANSPDB,ANSAMB
        CHARACTER*3 SECNM,SECNMN
        CHARACTER*4 ATNM,AN1,AN2,AN3,AN4,AC5,AO5,AC4,AC3,AO3,AP,AO4,
     1              AC2,AC1,AO1
        DOUBLE PRECISION DIST,DIS,EL
         EQUIVALENCE (XYZ(1,1),CR(1,1))
        DATA EE/'2*EXO-3*ENDO','3*ENDO','3*ENDO-4*EXO','4*EXO','4*EXO-O*
     1ENDO','O*ENDO','O*ENDO-1*EXO','1*EXO','1*EXO-2*ENDO','2*ENDO', '2*
     2ENDO-3*EXO','3*EXO','3*EXO-4*ENDO','4*ENDO','4*ENDO-O*EXO','O*EXO'
     3,'O*EXO-1*ENDO','1*ENDO','1*ENDO-2*EXO','2*EXO','2*EXO-3*ENDO'/

 1111   FORMAT(A20)
C
C         FINDING THE CONNECTIVITY MATRIX OF THE ATOMS
C
C        OPEN(UNIT=14,FILE='CONNECT.OUT',STATUS='NEW')
      DATA AP,AO5,AC5,AC4,AC3,AO3/'P   ','O5'' ','C5'' ','C4'' ','C3'' '
     1                           ,'O3'' '/
 	DATA AC2,AC1,AO1,AO4/'C2'' ','C1'' ','O1'' ','O4'' '/
	CONV = 3.14159/180.0
	DO KS=1,NSEG
	     ALPHA(KS) = 0.0
	     BETA(KS) = 0.0
	     GAMA(KS) = 0.0
	     DELTA(KS) = 0.0
	     EPSILN(KS) = 0.0
	     ZETA(KS) = 0.0
	     CHI(KS) = 0.0	
	     ANEU0(KS) = 0.0
	     ANEU1(KS) = 0.0
	     ANEU2(KS) = 0.0
	     ANEU3(KS) = 0.0
             ANEU4(KS) = 0.0
        END DO
	IF(ANSPDB.NE.'Y') THEN
	   WRITE(6,*) 'Is the file written in AMBER format? (Y/N)'
	   READ(5,'(A1)') ANSAMB
	   IF(ANSAMB.EQ.'Y'.OR.ANSAMB.EQ.'y') THEN
	      KSEG = 0
	      DO KRES=1,NSEG
	         IF(SECNM(KRES).EQ.'HB '.OR.SECNM(KRES).EQ.'POM') THEN
	            KSEG = KSEG + 1
	            SECNMN(KSEG) = SECNM(KRES+1)
 	            ISTSN(KSEG) = ISTS(KRES)
	            IENSN(KSEG) = IENS(KRES+1)
	         END IF
	         IF(SECNM(KRES).EQ.'HE ') THEN
	            SECNMN(KSEG) = SECNM(KRES-1)
	            IENSN(KSEG) = IENS(KRES)
	         END IF
	      END DO
	      DO KRES=1,NSEG
	         ISTS(KRES) = 0
	         IENS(KRES) = 0
	         SECNM(KRES) = '   '
	      END DO
	      DO KRES=1,KSEG
C
	       ISTS(KRES) = ISTSN(KRES)
	       IENS(KRES) = IENSN(KRES)
	       SECNM(KRES) = SECNMN(KRES)
	    END DO
	    NSEG = KSEG
	  END IF
	END IF
           DO I=1,NAT
              DO J=1,4
                 NN(I,J) = 0
              END DO
           END DO
C
        DO 10 I=1,NAT
            KJ=0
            DO J=1,(I-1)
	      DO KN=1,4
 	        IF(I.EQ.NN(J,KN)) THEN          
	           KJ = KJ + 1
	           NN(I,KJ) = J
	        END IF
              END DO
	    END DO
	    J = I + 1
	    DO WHILE(J.LE.NAT.AND.J.LE.(I+60))
                   DIST=0.0
                   DO K=1,3
                       DIST=DIST+(XYZ(I,K)-XYZ(J,K))*(XYZ(I,K)-XYZ(J,K))
                   END  DO
                   IF(DIST .LT. 2.73  .AND. KJ .LT. 4) THEN
                        KJ=KJ+1
                        NN(I,KJ)=J
                   END IF
	           J = J + 1
            END DO
  10    CONTINUE
C          DO  I=1,NAT
C            WRITE(14,999)I,ATNM(I),(NN(I,KJ),KJ=1,4)
C 999    FORMAT(3X,I3,2X,A4,4(2X,I3))
C          END DO
C
        DO 100 I=1,NAT
            AN1 = ATNM(I)
           DO 200 K=1,4
              N1=NN(I,K)
              IF(N1 .EQ. 0)GO TO 200
 	      DO KS=1,NSEG
	         IF(N1.GE.ISTS(KS).AND.N1.LT.IENS(KS)) THEN
  	            IUNIT = KS
                 END IF
              END DO
            AN2 = ATNM(N1)
              DO  300 L=1,4
                  N2=NN(N1,l)
                 IF(N2.LE.N1.OR.N2.EQ.I)GO TO 300
                    AN3 = ATNM(N2)
                  DO  400 M=1,4
                       N3=NN(N2,M)
                       IF(N3.EQ.0.OR.N3.EQ.N1)GO TO 400
                          AN4 = ATNM(N3)
                       DO IJ=1,3
                           XYZ1(IJ)=XYZ(I,IJ)
                           XYZ2(IJ)=XYZ(N1,IJ)
                           XYZ3(IJ)=XYZ(N2,IJ)
                           XYZ4(IJ)=XYZ(N3,IJ)
                       END DO
      IF(AN1.EQ.AO5.AND.AN2.EQ.AC5.AND.AN3.EQ.AC4.AND.AN4.EQ.AC3) THEN
              CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
              GAMA(IUNIT) = TAU
      END IF
      IF(AN1.EQ.AC3.AND.AN2.EQ.AC4.AND.AN3.EQ.AC5.AND.AN4.EQ.AO5) THEN
	      CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	      GAMA(IUNIT) = TAU
      END IF
      IF(AN1.EQ.AP.AND.AN2.EQ.AO5.AND.AN3.EQ.AC5.AND.AN4.EQ.AC4) THEN
              CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
              BETA(IUNIT) = TAU
      END IF
      IF(AN1.EQ.AC4.AND.AN2.EQ.AC5.AND.AN3.EQ.AO5.AND.AN4.EQ.AP) THEN
	      CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	      BETA(IUNIT) = TAU
      END IF
      IF(AN1.EQ.AC5.AND.AN2.EQ.AC4.AND.AN3.EQ.AC3.AND.AN4.EQ.AO3) THEN
              CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
              DELTA(IUNIT) = TAU
      END IF
      IF(AN1.EQ.AO3.AND.AN2.EQ.AC3.AND.AN3.EQ.AC4.AND.AN4.EQ.AC5) THEN
	      CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	      DELTA(IUNIT) = TAU
      END IF
      IF(AN1.EQ.AC4.AND.AN2.EQ.AC3.AND.AN3.EQ.AO3.AND.AN4.EQ.AP) THEN
              CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
              EPSILN(IUNIT) = TAU
      END IF
      IF(AN1.EQ.AP.AND.AN2.EQ.AO3.AND.AN3.EQ.AC3.AND.AN4.EQ.AC4) THEN
	      CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	      EPSILN(IUNIT) = TAU
      END IF
      IF(AN1.EQ.AC3.AND.AN2.EQ.AO3.AND.AN3.EQ.AP.AND.AN4.EQ.AO5) THEN
              CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	      ZETA(IUNIT) = TAU
      END IF
      IF(AN1.EQ.AO5.AND.AN2.EQ.AP.AND.AN3.EQ.AO3.AND.AN4.EQ.AC3) THEN
	      CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	      ZETA(IUNIT) = TAU
      END IF
      IF(AN1.EQ.AO3.AND.AN2.EQ.AP.AND.AN3.EQ.AO5.AND.AN4.EQ.AC5) THEN
              CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
              ALPHA(IUNIT) = TAU
      END IF
      IF(AN1.EQ.AC5.AND.AN2.EQ.AO5.AND.AN3.EQ.AP.AND.AN4.EQ.AO3) THEN
	      CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	      ALPHA(IUNIT) = TAU
      END IF
      IF((AN1.EQ.'O4'' '.OR.AN1.EQ.'O1'' ').AND.AN2.EQ.'C1'' ') THEN
             IF(AN3.EQ.'N9  '.AND.AN4.EQ.'C4  ') THEN
                 CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
                 CHI(IUNIT) = TAU 
             END IF
             IF(AN3.EQ.'N1  '.AND.AN4.EQ.'C2  ') THEN
                 CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
                 CHI(IUNIT) = TAU 
             END IF
      END IF
	IF(AN1.EQ.AC4.AND.AN2.EQ.AO1.AND.AN3.EQ.AC1.AND.AN4.EQ.AC2)THEN
              CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
    	      ANEU0(IUNIT) = TAU
	END IF
	IF(AN1.EQ.AC2.AND.AN2.EQ.AC1.AND.AN3.EQ.AO1.AND.AN4.EQ.AC4) THEN
	      CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	      ANEU0(IUNIT) = TAU
	END IF
	IF(AN1.EQ.AC4.AND.AN2.EQ.AO4.AND.AN3.EQ.AC1.AND.AN4.EQ.AC2)THEN
              CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	      ANEU0(IUNIT) = TAU
	END IF
	IF(AN1.EQ.AC2.AND.AN2.EQ.AC1.AND.AN3.EQ.AO4.AND.AN4.EQ.AC4) THEN
	      CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	      ANEU0(IUNIT) = TAU
	END IF
	IF(AN1.EQ.AC3.AND.AN2.EQ.AC2.AND.AN3.EQ.AC1.AND.AN4.EQ.AO1)THEN
              CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	      ANEU1(IUNIT) = TAU
	END IF
	IF(AN1.EQ.AO1.AND.AN2.EQ.AC1.AND.AN3.EQ.AC2.AND.AN4.EQ.AC3) THEN
	      CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	      ANEU1(IUNIT) = TAU
	END IF
	IF(AN1.EQ.AC3.AND.AN2.EQ.AC2.AND.AN3.EQ.AC1.AND.AN4.EQ.AO4)THEN
              CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	      ANEU1(IUNIT) = TAU
	END IF
	IF(AN1.EQ.AO1.AND.AN2.EQ.AC1.AND.AN3.EQ.AC2.AND.AN4.EQ.AC3)THEN
           CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	   ANEU1(IUNIT) = TAU
	END IF
	IF(AN1.EQ.AO4.AND.AN2.EQ.AC1.AND.AN3.EQ.AC2.AND.AN4.EQ.AC3)THEN
          CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	  ANEU1(IUNIT) = TAU
	END IF
	IF(AN1.EQ.AC4.AND.AN2.EQ.AC3.AND.AN3.EQ.AC2.AND.AN4.EQ.AC1)THEN
           CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	   ANEU2(IUNIT) = TAU
	END IF
	IF(AN1.EQ.AC1.AND.AN2.EQ.AC2.AND.AN3.EQ.AC3.AND.AN4.EQ.AC4) THEN
	      CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	      ANEU2(IUNIT) = TAU
	END IF
	IF(AN1.EQ.AO4.AND.AN2.EQ.AC4.AND.AN3.EQ.AC3.AND.AN4.EQ.AC2)THEN
           CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	   ANEU3(IUNIT) = TAU
	END IF
	IF(AN1.EQ.AC2.AND.AN2.EQ.AC3.AND.AN3.EQ.AC4.AND.AN4.EQ.AO4) THEN
	      CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	      ANEU3(IUNIT) = TAU
	END IF
	IF(AN1.EQ.AO1.AND.AN2.EQ.AC4.AND.AN3.EQ.AC3.AND.AN4.EQ.AC2)THEN
           CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	   ANEU3(IUNIT) = TAU
	END IF
	IF(AN1.EQ.AC2.AND.AN2.EQ.AC3.AND.AN3.EQ.AC4.AND.AN4.EQ.AO1) THEN
	      CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	      ANEU3(IUNIT) = TAU
	END IF
	IF(AN1.EQ.AC3.AND.AN2.EQ.AC4.AND.AN3.EQ.AO4.AND.AN4.EQ.AC1)THEN
           CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	   ANEU4(IUNIT) = TAU
	END IF
	IF(AN1.EQ.AC1.AND.AN2.EQ.AO4.AND.AN3.EQ.AC4.AND.AN4.EQ.AC3) THEN
	      CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	      ANEU4(IUNIT) = TAU
	END IF
	IF(AN1.EQ.AC3.AND.AN2.EQ.AC4.AND.AN3.EQ.AO1.AND.AN4.EQ.AC1)THEN
           CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	   ANEU4(IUNIT) = TAU
	END IF
	IF(AN1.EQ.AC1.AND.AN2.EQ.AO1.AND.AN3.EQ.AC4.AND.AN4.EQ.AC3) THEN
	      CALL TORAN(XYZ1,XYZ2,XYZ3,XYZ4,TAU)
	      ANEU4(IUNIT) = TAU
	END IF
C
 400    CONTINUE
 300    CONTINUE
 200    CONTINUE
 100    CONTINUE        
C
        WRITE(11,110)
 110    FORMAT('1',2X,'BACKBONE AND GLYCOSIDIC TORSION ANGLES'/) 
  	WRITE(11,111)
 111    FORMAT(7X,' P-O5''  O5''-C5''  C5''-C4''  C4''-C3''  C3''-O3'' 
     1 O3''- P    C1''-N')
        WRITE(11,222)
 222    FORMAT(7X,'ALPHA     BETA    GAMMA    DELTA     EPS      ZETA   
     1  CHI'/)
        DO KS=1,NSEG
	 WRITE(11,'(3X,7F9.1)') ALPHA(KS),BETA(KS),GAMA(KS),DELTA(KS),
     1                          EPSILN(KS),ZETA(KS),CHI(KS)
C
         IF(ALPHA(KS) .LT. 0.00)ALPHA(KS) = ALPHA(KS) + 360.0
         IF(BETA(KS) .LT. 0.00)BETA(KS) = BETA(KS) + 360.0
         IF(GAMA(KS) .LT. 0.00)GAMA(KS) = GAMA(KS) + 360.0
         IF(DELTA(KS) .LT. 0.00)DELTA(KS) = DELTA(KS) + 360.0
         IF(EPSILN(KS) .LT. 0.00)EPSILN(KS) = EPSILN(KS) + 360.0
         IF(ZETA(KS) .LT. 0.00)ZETA(KS) = ZETA(KS) + 360.0
         IF(CHI(KS) .LT. 0.00)CHI(KS) = CHI(KS) + 360.0
	END DO
       NSEG1 = NSEG - 1
      IF(ANSDB .EQ. 'Y')NSEG1 = NSEG - 2
        CALL STATIS(ALPHA,NSEG,AVALF,STDALF)
        CALL STATIS(BETA,NSEG,AVBET,STDBET)
        CALL STATIS(GAMA,NSEG,AVGAM,STDGAM)
        CALL STATIS(DELTA,NSEG,AVDEL,STDDEL)
        CALL STATIS(EPSILN,NSEG,AVEPS,STDEPS)
        CALL STATIS(ZETA,NSEG,AVZET,STDZET)
        CALL STATIS(CHI,NSEG,AVCHI,STDCHI)
        IF(AVALF .GT. 180.0)AVALF = AVALF - 360.0
        IF(AVBET .GT. 180.0)AVBET = AVBET - 360.0
        IF(AVGAM .GT. 180.0)AVGAM = AVGAM - 360.0
        IF(AVDEL .GT. 180.0)AVDEL = AVDEL - 360.0
        IF(AVEPS .GT. 180.0)AVEPS = AVEPS - 360.0
        IF(AVZET .GT. 180.0)AVZET = AVZET - 360.0
        IF(AVCHI .GT. 180.0)AVCHI = AVCHI - 360.0
       WRITE(11,120)AVALF,AVBET,AVGAM,AVDEL,AVEPS,AVZET,AVCHI
 120   FORMAT(/,1X,'Avg.',F7.1,6F9.1)
 121   FORMAT(1X,'S.D.',F7.1,6F9.1)
       WRITE(11,121)STDALF,STDBET,STDGAM,STDDEL,STDEPS,STDZET,STDCHI
C      
	WRITE(11,112)
 112    FORMAT(//5X,' O4''-C1''  C1''-C2''  C2''-C3''  C3''-C4''  C4''-O
     14''    AMPL    PHASE   CLASS'/)
C
	SS3672 = SIN(36.0 * CONV) + SIN(72.0 * CONV)
C
	 DO KS=1,NSEG
             TANPN = (ANEU4(KS) + ANEU1(KS)) - (ANEU3(KS) + ANEU0(KS))
	   IF(ANEU2(KS).NE.0.000E0) THEN
             TANPD = 2.0 * ANEU2(KS) * SS3672
	     PHASE(KS) = ATAN2(TANPN,TANPD) 
	     ANEUMAX(KS) = ANEU2(KS)/COS(PHASE(KS))
	     PHASE(KS) = PHASE(KS) / CONV
           ELSE
             PHASE(KS) = 00.000
 	   END IF
C  --  PHASE ANGLE IS IN DEGREES NOW
         PH = PHASE(KS)
         IF(PH .LT. 0.0)PH = 360.0 + PHASE(KS)
         CLASS = EE(INT((PH + 9)/18) + 1) 
         WRITE(11,'(3X,7F9.1,2X,A12)') ANEU0(KS),ANEU1(KS),
     1     ANEU2(KS),ANEU3(KS),ANEU4(KS),ANEUMAX(KS),PHASE(KS),CLASS
         END DO
C
       CALL STATIS(ANEU0,NSEG,AVN0,STDN0)
       CALL STATIS(ANEU1,NSEG,AVN1,STDN1)
       CALL STATIS(ANEU2,NSEG,AVN2,STDN2)
       CALL STATIS(ANEU3,NSEG,AVN3,STDN3)
       CALL STATIS(ANEU4,NSEG,AVN4,STDN4)
       CALL STATIS(ANEUMAX,NSEG,AVNMAX,STDNMAX)
       CALL STATIS(PHASE,NSEG,AVPHS,STDPHS)
         WRITE(11,120)AVN0,AVN1,AVN2,AVN3,AVN4,AVNMAX,AVPHS
         WRITE(11,121)STDN0,STDN1,STDN2,STDN3,STDN4,STDNMAX,STDPHS
C
	RETURN
        END
C
      SUBROUTINE TORAN(A,B,C,D,CHI)
      DIMENSION A(3),B(3),C(3),D(3)
      DIMENSION WX(3),WY(3),WZ(3),EL(3),EM(3),EN(3)
      DO 100 J = 1,3
      WX(J)=A(J)-B(J)
      WY(J)=C(J)-B(J)
      WZ(J)=D(J)-B(J)
100   CONTINUE
      CALL UNITV(WX(1),WX(2),WX(3),WY(1),WY(2),WY(3),EL(1),EM(1),EN(1))
      CALL UNITV(WZ(1),WZ(2),WZ(3),WY(1),WY(2),WY(3),EL(2),EM(2),EN(2))
      CHI=ACOSF(EL(1)*EL(2)+EM(1)*EM(2)+EN(1)*EN(2))
      CALL UNITV(EL(2),EM(2),EN(2),EL(1),EM(1),EN(1),EL(3),EM(3),EN(3))
      CHECK=EL(3)*WY(1)+EM(3)*WY(2)+EN(3)*WY(3)
      IF(CHECK.LT.0.0) RETURN
      CHI=0.0-CHI
      IF(CHI.GT.180.0) THEN
          CHI = CHI - 360.0
      END IF
      RETURN
      END

C***********************************************************************

      SUBROUTINE UNITV(A1,A2,A3,B1,B2,B3,EL,EM,EN)
      C1 = A2*B3-A3*B2
      C2 = A3*B1-A1*B3
      C3 = A1*B2-A2*B1
      FAC = 0.0
      FAC1 =SQRTF(C1*C1+C2*C2+C3*C3)
      IF(FAC1.LE.0.00001) GO TO 100
      FAC = 1.0/FAC1
100   CONTINUE
      EL =C1*FAC
      EM =C2*FAC
      EN =C3*FAC
      RETURN
      END
C
      FUNCTION ACOSF(ANG)
1     FORMAT('  Cos (Theta) is greater than 1 !',F12.8)
      IF(ABS(ANG).LE.1.0) GO TO 102
C     WRITE(1,1) ANG
      IF (ABS(ANG).GT.1.000001) GO TO 102
      ACOSF=0.0
      IF (ANG.LT.0.0) ACOSF = 180.0
      RETURN
102   CONTINUE
      ACOSF=ACOS(ANG)*180.0/3.14159265
      RETURN
      END
C *********************************************************************
C
      FUNCTION SQRTF(X)
1     FORMAT(2X,'Square root of a negative quantity!',E16.8)
      IF (X.LT.0.0) GO TO 101
      SQRTF = SQRT(X)
      RETURN
101   CONTINUE
C     WRITE(1,1) X
      SQRTF = 0.0
      RETURN
      END

C***********************************************************************
C
      SUBROUTINE LINEFIT(DCS,NDATA) 
      DOUBLE PRECISION DCS,DCT
      LOGICAL RLAMDA
      COMMON /DISPL/DXLOC(99),DYLOC(99),DZLOC(99),KKTI,
     1      BSCENTR(99,3),HORIG(99,3)
	COMMON /BASINF/BASE(99)
      DIMENSION DX(99),DY(99),DZ(99),AL(3),DIS(3),DCS(3),
     1          EL(3),EM(3),EN(3),DCT(3)
	COMMON /CENTER/CENTRS(99,3)
      CHARACTER*3 BASE
      CHARACTER*1 ANSO
C
      RLAMDA = .TRUE.
      XSTAR = 0.0
      YSTAR = 0.0
      ZSTAR = 0.0
      IF(NDATA.GT.2) THEN
        DO I=1,NDATA
          XSTAR = XSTAR + CENTRS(I,1)
          YSTAR = YSTAR + CENTRS(I,2)
          ZSTAR = ZSTAR + CENTRS(I,3)
        END DO
        XSTAR = XSTAR/NDATA
        YSTAR = YSTAR/NDATA
        ZSTAR = ZSTAR/NDATA
        AM11 = 0.0
        AM12 = 0.0
        AM13 = 0.0
        AM22 = 0.0
        AM23 = 0.0
        AM33 = 0.0
        DO I=1,NDATA
           DX(I) = CENTRS(I,1) - XSTAR
           DY(I) = CENTRS(I,2) - YSTAR
           DZ(I) = CENTRS(I,3) - ZSTAR
           AM11 = AM11 + DX(I) * DX(I)
           AM12 = AM12 + DX(I) * DY(I)
           AM13 = AM13 + DX(I) * DZ(I)
           AM22 = AM22 + DY(I) * DY(I)
           AM23 = AM23 + DY(I) * DZ(I)
           AM33 = AM33 + DZ(I) * DZ(I)
      END DO
      A = AM11
      B = AM12
      C = AM13
      D = AM22
      E = AM23
      F = AM33
      C1 = A*D*F + 2.0*B*C*E - A*E*E - B*B*F - C*C*D
      C2 = A*D + A*F + D*F - B*B - C*C - E*E
      C3 = A + D + F
      ALAMDA = 0.0
      ALH = 5.0
      DO 100 I=1,1000
      IF (ABS(ALH).GT.0.000001) THEN
      AFPR = (C2 + 2.0*ALAMDA*C3 - 3.0 * ALAMDA * ALAMDA)
      AF = C1 - ALAMDA*C2 + ALAMDA*ALAMDA*C3 - ALAMDA*ALAMDA*ALAMDA
      ALH = -AF/AFPR
      ALAMDA = ALAMDA + ALH
      ELSE
         GO TO 200
      END IF
 100  CONTINUE
 200  CONTINUE
C     WRITE(6,*) 'FIRST EIGEN-VALUE IS',ALAMDA
      A1 = 1.0
      A2 = ALAMDA - C3
      A3 = C1/ALAMDA
      SSQR = A2 * A2 - 4.0 * A3
      IF(SSQR.GT.0.0) THEN
         SSQR = SQRT(SSQR)
         ALAMD2 = (-A2 + SSQR) * 0.5
         ALAMD3 = (-A2 - SSQR) * 0.5
C        WRITE(6,*) 'SECOND EIGEN-VALUE IS',ALAMD2
C        WRITE(6,*) 'THIRD EIGEN-VALUE IS',ALAMD3
      ELSE
         SSQR = SQRT(-SSQR)
	 RLAMDA = .FALSE.
         A2 = -A2
C        WRITE(6,*) 'SECOND EIGEN-VALUE IS',A2,'i',SSQR
         SSQR = -SSQR
C        WRITE(6,*) 'THIRD EIGEN-VALUE IS',A2,'i',SSQR
      END IF
      AL(1) = ALAMDA
      AL(2) = ALAMD2
      AL(3) = ALAMD3
      KTERM = 1
      IF(RLAMDA) KTERM = 3
      DO KT=1,KTERM
         ADIS = 0.0
         A1 = AM11 - AL(KT)
         B1 = AM12
         C1 = AM13
         A2 = B1
         B2 = AM22 - AL(KT)
         C2 = AM23
C
         ELT = (B1 * C2 - C1 * B2)
         EMT = (C1 * A2 - A1 * C2)
         ENT = (A1 * B2 - B1 * A2)
         CONS = SQRT(ELT * ELT + EMT * EMT + ENT * ENT)
         ELT = ELT/CONS
         EMT = EMT/CONS
         ENT = ENT/CONS
C     WRITE(6,*) 'ELT,EMT,ENT',ELT,EMT,ENT,CONS
         DO K=1,NDATA
           ADIS = ADIS + DX(K)*DX(K) + DY(K)* DY(K) + DZ(K)*DZ(K)
           ADIS = ADIS - (ELT*DX(K)+EMT*DY(K)+ENT*DZ(K))**2
         END DO
C     WRITE(6,*) 'MEAN DIST.',ADIS
         ADIS = SQRT(ADIS)
         DIS(KT) = ADIS
         EL(KT) = ELT
         EM(KT) = EMT
         EN(KT) = ENT
         IF(KT.GT.1) THEN
            IF(DIS(KT).LT.DIS(KT-1)) NN=KT
         END IF
       END DO
       IF(RLAMDA) THEN
        DCS(1) = EL(NN) * 1.0
        DCS(2) = EM(NN) * 1.0
        DCS(3) = EN(NN) * 1.0
       ELSE
        DCS(1) = EL(1) * 1.0
	DCS(2) = EM(1) * 1.0
	DCS(3) = EN(1) * 1.0
       END IF
      ELSE
	  DO KKN=1,3
	    DCT(KKN) = CENTRS(1,KKN) - CENTRS(2,KKN)
          END DO
	  CALL NORMAL(DCT(1),DCT(2),DCT(3),DCS(1),DCS(2),DCS(3))
      END IF
      IF(CENTRS(1,3).GT.CENTRS(NDATA,3)) THEN
CWRITE(6,*) 'MOLECULE GOES DOWN IN Z'
     	DCS(1) = -DCS(1)
	DCS(2) = -DCS(2)
	DCS(3) = -DCS(3)
      END IF
      RETURN
      END
C***********************************************************************
      FUNCTION ACOSD(X)
      IF(ABS(X).LT.1.0) THEN
         CONV=180.0/3.141592654
         ACOSD = ACOS(X) * CONV
      ELSE
C        WRITE(1,*) 'VALUE OF COS(THETA) > 1.0'
         ACOSD = 0.0
      END IF
      RETURN
      END
C
C-----------------------------------------------------------------------
 	SUBROUTINE RTHEPH(ANSDB)
         COMMON/ATINF/CR(2400,3),SECSQ(300),
     1             ISTS(300),IENS(300),NSEG,NAT
	COMMON /ATINFC/ATNM(2400),SECNM(300)
	DIMENSION RAD(99),PHI(99),ZH(99),XH(99),YH(99),PPD1(99),PPD2(99)
	CHARACTER*4 ATNM, ANSDB*1,ANS*1
	CHARACTER*3 SECNM
        CONV=180.0/3.141592654
C
	I = 0
	DO II=1,NAT
	   IF(ATNM(II)(1:1).EQ.'P') THEN
              I = I + 1
	      RAD(I) = SQRT(CR(II,1)*CR(II,1)+CR(II,2)*CR(II,2))
	      PHI(I) = ATAN2(CR(II,2),CR(II,1)) * CONV
	      ZH(I) = CR(II,3)
              XH(I) = CR(II,1)
              YH(I) = CR(II,2)
	   END IF
	END DO
	WRITE(11,1)
C
C    DELPHI, DELZ AND P--P DISTANCE CALC.

	KK=I

	 IF ((ANSDB .EQ. 'Y')) KK = I/2

         I1 = 1
         I2 = KK - 1

	 JJJ=1
	 IF ((ANSDB .EQ. 'Y')) JJJ=2

       DO JJ = 1,JJJ
        WRITE(11,99)JJ
          IP=0
         DO NP = I1,I2
          IP = IP + 1
          NP1 = NP + 1
           DELP = PHI(NP1) - PHI(NP)
           IF(DELP .GT. 180.0)DELP = DELP - 360.0
           IF(DELP .LT. -180.0 .AND. DELP .GT. -360.0)DELP=360.0+DELP
C
           DELZ = ZH(NP1)  - ZH(NP)
           PPD1(NP)=SQRT((XH(NP)-XH(NP1))**2 + (YH(NP)-YH(NP1))**2 +
     1     (ZH(NP)-ZH(NP1))**2)
          WRITE(11,100)IP,RAD(NP),PHI(NP),ZH(NP),DELP,DELZ,PPD1(NP)
         END DO
          WRITE(11,100)IP+1,RAD(NP1),PHI(NP1),ZH(NP1)
         I1 = KK +1
         I2 = I - 1
        END DO  

	 IF ((ANSDB .EQ. 'Y')) THEN
C    INTERCHAIN P--P DISTS.

        WRITE(11,101)
        DO NP = 1,KK
         KI = KK + 1
         II=0
          DO NP2 = KI,I
           II = II + 1
           PPD2(II)=SQRT((XH(NP)-XH(NP2))**2 + (YH(NP)-YH(NP2))**2 + 
     1     (ZH(NP)-ZH(NP2))**2)
          END DO
         WRITE(11,102)NP,(PPD2(IP),IP=1,II)
        END DO
	ENDIF

  1	FORMAT('1','PHOSPHATE ATOMIC COORDINATES & P--P DISTANCES'//
     111X,'R',6X, ' PHI',8X,'Z',8X,'DELP',6X,'DELZ',8X,'P--P',/)
  99   FORMAT(3X,'Strand :',I2)
 100   FORMAT(2X,'P',I2,'.',3F9.2,2X,2F9.2,3X,F9.2)
 101   FORMAT(//2X,' INTER-CHAIN P--P DISTANCES'/'   STRAND 2 -->'/' STR
     1AND 1')
 102   FORMAT(2X,'P',I2,'.',10F7.2,2X,/,6X,10F7.2)
 199   FORMAT( '1')
	RETURN
	END
C
	SUBROUTINE FRACTION(A,B,C,ALPHA,BETA,GAMA)
        COMMON/ATINF/CR(2400,3),SECSQ(300),
     1             ISTS(300),IENS(300),NSEG,NAT
	COMMON /ATINFC/ATNM(2400),SECNM(300)
	DIMENSION AMAT(3,3),TEM(3)
	CHARACTER*3 SECNM
	CHARACTER*4 ATNM
        CONV=180.0/3.141592654
	ALPHA = ALPHA/CONV
	BETA = BETA/CONV	
	GAMA = GAMA/CONV
	AMAT(1,1)=1.0
	AMAT(1,2) = COS(GAMA)
	AMAT(1,3) = COS(BETA)
	AMAT(2,1) = 0.0
	AMAT(2,2) = SIN(GAMA)
	P1 = (COS(ALPHA) - COS(BETA) * COS(GAMA))/SIN(GAMA)	
	P2 = SQRT(SIN(BETA) * SIN(BETA) - P1 * P1)
	AMAT(2,3) = P1
	AMAT(3,1) = 0.0
	AMAT(3,2) = 0.0	
	AMAT(3,3) = P2
	DO KK=1,NAT
	   CR(KK,1) = CR(KK,1) * A
 	   CR(KK,2) = CR(KK,2) * B
  	   CR(KK,3) = CR(KK,3) * C
	   CALL MATMUL(AMAT,CR(KK,1),CR(KK,2),CR(KK,3),TEM)
	   DO KKK=1,3
	      CR(KK,KKK) = TEM(KKK)
	   END DO
	END DO
	RETURN
	END
C
	SUBROUTINE REORNT(IBAPR,NB,ANSORN,NPTS)
	DIMENSION IBAPR(200),COOR1(3),COOR2(3)
	COMMON /CENTER/CENTRS(99,3)
      COMMON/ATINF/CR(2400,3),SECSQ(300),ISTS(300)
     1            ,IENS(300),NSEG,NAT
	COMMON /ATINFC/ATNM(2400),SECNM(300)
      COMMON /DISPL/DXLOC(99),DYLOC(99),DZLOC(99),KKTI,
     1      BSCNTR(99,3),HORIG(99,3)
	COMMON /BASINF/BASE(99)
	CHARACTER*3 SECNM,BASE
	CHARACTER*4 ATNM,ATOM,ANSORN
	CHARACTER*1 ANSBC,ANSBB
 1	CONTINUE
	IF(ANSORN.EQ.'    ') THEN
	   WRITE(6,91)
 	   WRITE(8,91)
  91	FORMAT(/' Should the GLOBAL AXIS be fitted to:'/' 1. Local HELIX
     1 ORIGINS (Type O)'/' 2. Base-Pair CENTERS (Type C)'/' 3. Any Backb
     2one atom (Type the atom name)')
 	   READ(5,'(A4)') ATOM
	   WRITE(8,'(5X,A4)') ATOM
	   ANSORN = ATOM
	END IF
C
	IF(ANSORN.EQ.'C   ') THEN
	   DO KK=1,NB + 1
	       DO KOO=1,3
	          CENTRS(KK,KOO) = BSCNTR(KK,KOO)
	       END DO
	   END DO
	   NPTS = NB + 1
	ELSE IF(ANSORN.EQ.'O   ') THEN
           DO KS =1,NB
              DO KOO=1,3
                 CENTRS(KS,KOO)=HORIG(KS,KOO)
              END DO
           END DO
	   NPTS = NB
	ELSE 
	  ATOM = ANSORN
	  ICNT = 0
          DO IS=1,NB
	     IFIRST = IBAPR(IS*2-1)
	     ISECND = IBAPR(IS*2)
             NDOUB1 = 0
	     DO IATM=ISTS(IFIRST),IENS(IFIRST)
	        IF(ATOM.EQ.ATNM(IATM)) THEN
	           ICNT = ICNT + 1
		   DO KOO=1,3
	              CENTRS(ICNT,KOO) = CR(IATM,KOO)
	           END DO
	        END IF
      	     END DO
	     DO IATM=ISTS(ISECND),IENS(ISECND)
	        IF(ATOM.EQ.ATNM(IATM)) THEN
	           ICNT  = ICNT + 1
	           DO KOO=1,3
		      CENTRS(ICNT,KOO) = CR(IATM,KOO)
	 	   END DO
	        END IF
	     END DO
 	   END DO	   
  	   ILAST = IBAPR((NB+1)*2-1)
	   ILST = IBAPR((NB+1)*2)
	   IF(IENS(ILST).LE.NAT) THEN
	      DO IATM=ISTS(ILAST),IENS(ILAST)
	        IF(ATOM.EQ.ATNM(IATM)) THEN
	           ICNT = ICNT + 1
	           DO KOO=1,3
	              CENTRS(ICNT,KOO) = CR(IATM,KOO)
	           END DO
	        END IF
      	     END DO
	     DO IATM=ISTS(ILST),IENS(ILST)
	        IF(ATOM.EQ.ATNM(IATM)) THEN
	           ICNT = ICNT + 1
	           DO KOO=1,3
		      CENTRS(ICNT,KOO) = CR(IATM,KOO)
	 	   END DO
	        END IF
	     END DO
 	  END IF
	     NPTS = ICNT
 	END IF

	RETURN
	END
  
	SUBROUTINE LINETM(DC,NDATA)
	DIMENSION VECTRS(3,99),DVA(99),NXA(99),DC(3),VZ(3)
	COMMON /CENTER/CENTRS(99,3)
	DOUBLE PRECISION DC,ELA,EMA,ENA
C
	DO KN=2,NDATA,2
	   NXA(KN) = KN
	   NXA(KN-1) = KN - 1
	   DO KOO=1,3
	     VECTRS(KOO,KN-1) = CENTRS(KN+1,KOO) - CENTRS(KN-1,KOO)
	     VECTRS(KOO,KN) = -CENTRS(KN,KOO) + CENTRS(KN+2,KOO)
	   END DO
	END DO
	CALL PLANE(VECTRS,(NDATA-2),NXA,0,DVA,ELA,EMA,ENA,DET,DVSX)
	DC(1) = ELA
	DC(2) = EMA
	DC(3) = ENA
	RETURN
	END
C
	SUBROUTINE POINSN(NDATA,DX,DY)
	COMMON /CENTER/CENTRS(99,3)
	DX = 0.0
	DY = 0.0
	DO I=1,(NDATA-4)
	   AM1 = (CENTRS(I,2)-CENTRS(I+1,2))/(CENTRS(I,1)-CENTRS(I+1,1))
	   AM1 = -1.0/AM1
	   POI11 = (CENTRS(I,1) + CENTRS(I+1,1)) * 0.5
	   POI21 = (CENTRS(I,2) + CENTRS(I+1,2)) * 0.5
	   AM2=(CENTRS(I+1,2)-CENTRS(I+2,2))/(CENTRS(I+1,1)-CENTRS(I+2,
     11))
	   AM2 = -1.0/AM2
	   POI12 = (CENTRS(I+1,1) + CENTRS(I+2,1)) * 0.5
	   POI22 = (CENTRS(I+1,2) + CENTRS(I+2,2)) * 0.5
	   C1 = POI21 - AM1 * POI11
	   C2 = POI22 - AM2 * POI12
	   POINT1 = (C1 - C2)/(AM2 - AM1)
	   POINT2 = (AM2 * C1 - AM1 * C2)/(AM2 - AM1)
	   DX = DX + POINT1
	   DY = DY + POINT2
	END DO
	DX = DX/(NDATA - 4)
	DY = DY/(NDATA - 4)
	RETURN
	END
C
	SUBROUTINE POINDB(NDATA,DX,DY)
	COMMON /CENTER/CENTRS(99,3)
C
	DX = 0.0
	DY = 0.0
	NVEC = NDATA/2
	DO I=1,(NDATA-2),2
         AM1=(CENTRS(I,2)-CENTRS(I+1,2))/(CENTRS(I,1)-CENTRS(I+1,1))
	 AM1 = -1.0/AM1
	 POI11 = (CENTRS(I,1) + CENTRS(I+1,1)) * 0.5
	 POI21 = (CENTRS(I,2) + CENTRS(I+1,2)) * 0.5
	 AM2=(CENTRS(I+2,2)-CENTRS(I+3,2))/(CENTRS(I+2,1)-CENTRS(I+3,1))
	 AM2 = -1.0/AM2
	 POI12 = (CENTRS(I+2,1) + CENTRS(I+3,1)) * 0.5
	 POI22 = (CENTRS(I+2,2) + CENTRS(I+3,2)) * 0.5
	 C1 = POI21 - AM1 * POI11
	 C2 = POI22 - AM2 * POI12
	 POINT1 = (C1 - C2)/(AM2 - AM1)
	 POINT2 = (AM2 * C1 - AM1 * C2)/(AM2 - AM1)
	 DX = DX + POINT1
	 DY = DY + POINT2
	END DO
	DX = DX/(NVEC - 1)
	DY = DY/(NVEC - 1)
	RETURN
	END
	SUBROUTINE LINESN(DC,NDATA)
	DIMENSION VECTRS(3,99),DVA(99),NXA(99),DC(3),VZ(3)
	COMMON /CENTER/CENTRS(99,3)
	DOUBLE PRECISION DC,ELA,EMA,ENA
C
	DO KN=2,NDATA
	   NXA(KN-1) = KN - 1
	   DO KOO=1,3
	     VECTRS(KOO,KN-1) = CENTRS(KN,KOO) - CENTRS(KN-1,KOO)
	   END DO
	END DO
	CALL PLANE(VECTRS,(NDATA-1),NXA,0,DVA,ELA,EMA,ENA,DET,DVSY)
	DC(1) = ELA
	DC(2) = EMA
	DC(3) = ENA
	RETURN
	END