src/contrib/total.f

   1 C Chemical shift calculation, to read in multiple NMR structures
   2 C (with protons) and calculate sd for each
   3 C Will also read Xray structures and add protons
   4 C Current limit is 5000 heavy atoms 3000 protons and 50 rings
   5 C Author - Mike Williamson Jan 94
   6 C see M P Williamson and T Asakura, J Magn Reson Ser B 101 63-71 1993
   7         DIMENSION SHIFT(3000)   !Shift values for each proton
   8         DIMENSION ANISCO(3000),ANISCN(3000),SHIFTE(3000)
   9         DIMENSION sum(3000),exptln(3000)
  10         DIMENSION EXPTL(3000)   !Exptl shifts (external file)
  11         DIMENSION IRES(5000),RES(5000),ANAME(5000),FINAL(3000)
  12         DIMENSION X(5000), Y(5000), Z(5000) !Heavy atom input
  13         DIMENSION IRESH(3000),RESH(3000),ANAMEH(3000)
  14         DIMENSION XH(3000), YH(3000), ZH(3000) !H atom input
  15         DIMENSION vx(3),vy(3),vz(3),CEN(3),rh(3),pregam(3),yy(3)
  16 C (Anisotropy tensor plus anisotropy centre)
  17         DIMENSION vca(3),vc(3),vo(3),vn(3),calctemp(3000)
  18         DIMENSION datres(179),datnam(179),datshift(179)
  19         DIMENSION RCNET(3000)
  20         DIMENSION IACODE(50),IRATPT(9,50),IRATS(50)
  21         DIMENSION shiftt(3000),exptlt(3000),VCHN(3)
  22         dimension rn(3,50),cent(3,50),signr(50)
  23         dimension iexp(2000),atexp(2000),pexp(2000),eexp(2000)
  24         CHARACTER ANAME*4,RES*3,FN1*60,FN3*40,junk*80
  25         CHARACTER ans*1,FN4*40,pexp*3,resh*3,anameh*4,prott*3
  26         CHARACTER datres*3,datnam*4,YN*1,atexp*4,FN9*60,RFILE*60
  27         integer AT1,AT2,AT3,INDX(3000)
  28         COMMON x,y,z,xh,yh,zh
  29 C Set values for anisotropies, atomic charges and multiplication
  30 C factor for electric field: shift=-Ez*sigmaE
  31         DATA XCO1/-13.0/,XCO2/-4.0/,XCN1/-11.0/,XCN2/1.40/
  32         DATA sigmaE/0.60/,Qc/1.60/,Qn/-1.70/,Qo/-2.30/
  33         DATA Qhn/0.70/
  34 C NB The atomic charges are in esu
  35 C
  36 C******************FILE INPUT************************************
  37         iwarning=1
  38         type *,'Input file?'
  39         read(5,999)FN1
  40 997     format (I)
  41 98      FORMAT (I5)
  42 999     format(A)
  43         type *,'Output file?'
  44         read(5,999)FN3
  45         type *,'Calculate for HA[1], HN[2] or all protons[3]?'
  46         read(5,997)icalculate
  47         print *,'icalculate=',icalculate
  48 C NB Actually calculates all protons, just prints differently
  49 99      FORMAT(A80)
  50         type *,' Are you comparing calc. to experimental shifts? [N]'
  51         read(5,970)YN
  52 970     format(A)
  53         OPEN(1,file=FN1,STATUS='OLD',readonly)
  54         OPEN(3,file=FN3)
  55         type *,' Random file ?'
  56         read(5,999)RFILE
  57         open(4,file=RFILE,status='old',readonly)
  58         do 25 I=1,179           !Random coil shifts
  59         read(4,998) datres(I),datnam(I),datshift(I)
  60 998     format(A3,1X,A4,F5.2)
  61 25      continue
  62         if(yn.eq.'Y'.or.yn.eq.'y') then
  63         type *,' File containing experimental shifts?'
  64         read(5,999)FN4
  65         OPEN(UNIT=8,FILE=FN4,STATUS='OLD')
  66 C This file should contain 2-line header, no. of protons (I5) plus list
  67         type *,'Name of protein in this file?'
  68         read(5,992) prott
  69 992     format(A3)
  70         READ (8,99) JUNK
  71         READ (8,99) JUNK
  72         READ (8,98) Iprot
  73         do 900 II=1,Iprot
  74         read (8,996) iexp(ii),atexp(ii),pexp(ii),eexp(ii)
  75 996     format(I3,7X,A4,5X,A3,27X,F9.5)
  76 900     continue
  77         endif
  78         imodel=0
  79 776     iheavyatom=0
  80         iproton=0
  81 C PDB file should end with TER or END or preferably both
  82         do 10 I=1,100000
  83         READ (1,99,end=777) JUNK
  84         if (junk(1:5).eq.'MODEL') then
  85          read(junk(6:14),'(I9)') nmodel
  86          imodel=1
  87         endif
  88         if (junk(1:6).eq.'ENDMDL') goto 11
  89         if (junk(1:4).eq.'END') then
  90           if (imodel.eq.1) goto 777
  91           goto 11
  92         endif
  93 C Replace D by H for neutron structures
  94         if (junk(14:14).eq.'D') junk(14:14)='H'
  95         if (junk(1:4).eq.'TER ') goto 11
  96         if (junk(1:4).eq.'ATOM') then
  97          if(junk(27:27).ne.' ') print 601,junk(23:26)
  98 601       format('WARNING: substituted residue present at ',A4)
  99          if(junk(22:22).ne.' '.and.iwarning.eq.1) then
 100           iwarning=2
 101           type *,'WARNING: more than one chain present'
 102          endif
 103          if(junk(17:17).ne.' ') print 602,junk(23:26)
 104 602       format('WARNING: alternate conformations at residue ',A4)
 105 C Next line reads HN in as heavy atom, for Electric field calc.
 106          if((junk(14:14).ne.'H').or.(junk(13:15).eq.' H ')) then
 107          iheavyatom=iheavyatom+1
 108          read(junk(13:16),'(A4)') aname(iheavyatom)
 109          read(junk(18:20),'(A3)') res(iheavyatom)
 110          read(junk(23:26),'(I4)') ires(iheavyatom)
 111          read(junk(31:54),'(3F8.0)') x(iheavyatom),
 112      .     y(iheavyatom),z(iheavyatom)
 113          endif
 114          if(junk(14:14).eq.'H') then
 115          if((icalculate.eq.1).and.(junk(14:15).ne.'HA'))goto 10
 116          if((icalculate.eq.2).and.(junk(14:15).ne.'H '))goto 10
 117          iproton=iproton+1
 118          read(junk(13:16),'(A4)') anameh(iproton)
 119          read(junk(18:20),'(A3)') resh(iproton)
 120          read(junk(23:26),'(I4)') iresh(iproton)
 121          read(junk(31:38),'(F8.3)') xh(iproton)
 122          read(junk(39:46),'(F8.3)') yh(iproton)
 123          read(junk(47:54),'(F8.3)') zh(iproton)
 124          endif
 125         endif
 126 10       CONTINUE
 127 11       continue
 128 C To avoid calculating for water molecules
 129          if(res(1).eq.'WAT'.or.res(2).eq.'WAT') goto 777
 130 C Now see if protons are present and add them if not
 131 C
 132 C       Next Line hacked (DvdS, 12/94)
 133         if ((iproton.eq.0) .or. (icalculate.eq.3)) then
 134         type *,'Adding protons'
 135         type *,'Print out file with protons?'
 136         read(5,607)ANS
 137         call addprot(x,y,z,xh,yh,zh,Iheavyatom,Iproton,
 138      &   aname,anameh,res,resh,ires,iresh)
 139         if(ans.eq.'Y'.or.ans.eq.'y') then
 140          do 670 I=1,60
 141           if(FN1(I:I).eq.' ')goto 671
 142 670      continue
 143 671      ilength=I-1
 144           if(ilength.lt.40) then
 145            FN9=FN1(1:ilength-4)//'_protonated.pdb'
 146           else
 147            FN9=FN1(ilength-8:ilength-4)//'_protonated.pdn'
 148           endif
 149           print 669,FN9
 150 669       format('Output going to ',A)
 151           open(9,file=FN9)
 152           write(9,660)FN1
 153 660       format('REMARK  Generated from ',60A)
 154           iresstart=ires(1)
 155           iresend=ires(iheavyatom)
 156           iline=0
 157           do 668 icount=iresstart,iresend
 158           do 661 I=1,iheavyatom
 159           if(aname(I).eq.' H  ') goto 661
 160           if(ires(I).eq.icount) then
 161           iline=iline+1
 162           write(9,662)iline,aname(I),res(I),ires(I),x(I),y(I),z(I)
 163           endif
 164 662       format('ATOM',I7,1X,A4,1X,A3,2X,I4,4X,3F8.3)
 165 661       continue
 166           do 663 I=1,iproton
 167           if(iresh(I).eq.icount) then
 168           iline=iline+1
 169           write(9,662)iline,anameh(I),resh(I),iresh(I),xh(I),yh(I),zh(I)
 170           endif
 171 663       continue
 172 668       continue
 173           write(9,665)
 174           write(9,666)
 175 665       format('TER')
 176 666       format('END')
 177         end if
 178        end if
 179          do 20 I=1,iproton
 180          shift(I)=0.0           !initialise
 181          anisco(I)=0.0
 182          aniscn(I)=0.0
 183          shiftE(I)=0.0
 184 20       continue
 185 607     FORMAT(A1)
 186         type *,' All atoms read...initialising'
 187         if(imodel.eq.1) type 774, nmodel
 188 774     format(' Calculating shifts for model',I9)
 189 C***********Calculate number of aromatic residues, rearrange order of
 190 C ring atoms, and set pointers to line numbers of aromatic atoms
 191 C (NB each Trp counts as two rings)
 192         narom=0
 193         do 105 I=1,iheavyatom
 194         if((res(I).eq.'TRP'.or.res(I).eq.'TYR'.or.res(I).eq.'PHE'.
 195      1   or.res(I).eq.'HIS').and.(aname(I).eq.' CB ')) THEN
 196          narom = narom + 1
 197          do 102 Iring=1,6
 198           IRATPT(Iring,narom)=0
 199 102      continue
 200          IACODE(narom)=IRES(I)
 201           IF ((res(I).eq.'PHE').or.(res(I).eq.'TYR')) THEN
 202             IRATS(narom)=6
 203             do 101 K=1,20
 204             IF(ires(I+K).eq.ires(I).and.
 205      &       aname(I+K).eq.' CG ') IRATPT(1,narom)=I+K
 206             IF(ires(I+K).eq.ires(I).and.
 207      &       aname(I+K).eq.' CD1') IRATPT(2,narom)=I+K
 208             IF(ires(I+K).eq.ires(I).and.
 209      &       aname(I+K).eq.' CE1') IRATPT(3,narom)=I+K
 210             IF(ires(I+K).eq.ires(I).and.
 211      &       aname(I+K).eq.' CZ ') IRATPT(4,narom)=I+K
 212             IF(ires(I+K).eq.ires(I).and.
 213      &       aname(I+K).eq.' CE2') IRATPT(5,narom)=I+K
 214             IF(ires(I+K).eq.ires(I).and.
 215      &       aname(I+K).eq.' CD2') IRATPT(6,narom)=I+K
 216 101         continue
 217           ELSE IF (res(I).eq.'HIS') THEN
 218             IRATS(narom)=5
 219             do 103 K=1,20
 220             IF(ires(I+K).eq.ires(I).and.
 221      &       aname(I+K).eq.' CG ') IRATPT(1,narom)=I+K
 222             IF(ires(I+K).eq.ires(I).and.
 223      &       aname(I+K).eq.' CD2') IRATPT(2,narom)=I+K
 224             IF(ires(I+K).eq.ires(I).and.
 225      &       aname(I+K).eq.' NE2') IRATPT(3,narom)=I+K
 226             IF(ires(I+K).eq.ires(I).and.
 227      &       aname(I+K).eq.' CE1') IRATPT(4,narom)=I+K
 228             IF(ires(I+K).eq.ires(I).and.
 229      &       aname(I+K).eq.' ND1') IRATPT(5,narom)=I+K
 230 103         continue
 231           ELSE
 232 C Trp counts as two aromatic rings (5 then 6)
 233             IRATS(narom)=5
 234             do 104 K=1,20
 235             IF(ires(I+K).eq.ires(I).and.
 236      &       aname(I+K).eq.' CG ') IRATPT(1,narom)=I+K
 237             IF(ires(I+K).eq.ires(I).and.
 238      &       aname(I+K).eq.' CD1') IRATPT(2,narom)=I+K
 239             IF(ires(I+K).eq.ires(I).and.
 240      &       aname(I+K).eq.' NE1') IRATPT(3,narom)=I+K
 241             IF(ires(I+K).eq.ires(I).and.
 242      &       aname(I+K).eq.' CE2') IRATPT(4,narom)=I+K
 243             IF(ires(I+K).eq.ires(I).and.
 244      &       aname(I+K).eq.' CD2') IRATPT(5,narom)=I+K
 245 104         continue
 246             narom = narom + 1
 247             IACODE(narom)=IRES(I)
 248             IRATS(narom)=6
 249             do 106 K=1,25
 250             IF(ires(I+K).eq.ires(I).and.
 251      &       aname(I+K).eq.' CE2') IRATPT(1,narom)=I+K
 252             IF(ires(I+K).eq.ires(I).and.
 253      &       aname(I+K).eq.' CD2') IRATPT(2,narom)=I+K
 254             IF(ires(I+K).eq.ires(I).and.
 255      &       aname(I+K).eq.' CE3') IRATPT(3,narom)=I+K
 256             IF(ires(I+K).eq.ires(I).and.
 257      &       aname(I+K).eq.' CZ3') IRATPT(4,narom)=I+K
 258             IF(ires(I+K).eq.ires(I).and.
 259      &       aname(I+K).eq.' CH2') IRATPT(5,narom)=I+K
 260             IF(ires(I+K).eq.ires(I).and.
 261      &       aname(I+K).eq.' CZ2') IRATPT(6,narom)=I+K
 262 106         continue
 263           ENDIF
 264         endif   !End of loop for each aromatic residue
 265 105     continue
 266         do 107 J=1,narom
 267         if(iratpt(1,j).eq.0.or.iratpt(2,j).eq.0.or.iratpt(3,j).eq.0.
 268      &   or.iratpt(4,j).eq.0.or.iratpt(5,j).eq.0) print 960,
 269      &   ires(iacode(j)),res(ires(iacode(j)))
 270 107     continue
 271 960     format(1X,'Ring atom missing for ',I4,1X,A3)
 272         DO 115 L=1,iproton      !Initialise all ring current shifts to zero
 273           RCNET(L)=0.
 274 115     CONTINUE
 275 C ******************************************************************
 276         type *,' Calculating ring current shifts....'
 277 C
 278       DO 116 J=1,NAROM          !FOR EACH AROMATIC RESIDUE
 279 C Now set up rn (ring normal) and cent for each ring.
 280 C Determined using atoms 1, 3 and 5 only of ring.
 281 C Much of this is lifted from a version of AMBER provided by Dave Case
 282         call plane(IRATPT(1,j),IRATPT(3,j),IRATPT(5,j),x,y,z,rn(1,j),
 283      .               cent(1,j))
 284 c  From pts 1, 3 and 5 calculates rn(ring normal), drn (deriv. of
 285 c  ring normal) and centre of ring
 286 c
 287 c  --   check on signr of normal vector
 288 c
 289             signr(j) = 1.0
 290             d1cx = x(IRATPT(1,j)) - cent(1,j)
 291             d1cy = y(IRATPT(1,j)) - cent(2,j)
 292             d1cz = z(IRATPT(1,j)) - cent(3,j)
 293             d2cx = x(IRATPT(3,j)) - cent(1,j)
 294             d2cy = y(IRATPT(3,j)) - cent(2,j)
 295             d2cz = z(IRATPT(3,j)) - cent(3,j)
 296             vp1 = d1cy*d2cz - d1cz*d2cy
 297             vp2 = d1cz*d2cx - d1cx*d2cz
 298             vp3 = d1cx*d2cy - d1cy*d2cx
 299             if ((vp3*rn(3,j)+vp2*rn(2,j)+vp1*rn(1,j))
 300      .               .gt.0.0) signr(j) = -1.0
 301   119     DO 120 L=1,iproton    !For each proton
 302         ip=L
 303 C Next line includes self aromatic shift for HA only
 304          IF(IRESH(L).EQ.IRES(IRATPT(1,J)).AND.(ANAMEH(L)(2:3).NE.
 305      &     'HA'.AND.ANAMEH(L)(2:3).NE.'H ')) GOTO 120
 306 c
 307 c  --     skip rings whose centre is more than 15A away
 308 c
 309             relc = (xh(ip)-cent(1,j))**2 + (yh(ip)-cent(2,j))**2
 310      .            +(zh(ip)-cent(3,j))**2
 311             if (relc.gt.225.) go to 120
 312 c
 313 c  --   loop over pairs of bonded atoms in the ring
 314 c
 315           do 80 k=1,irats(j)
 316             kp1 = k + 1
 317             if (kp1.gt.irats(j)) kp1 = 1
 318 c
 319             call hm(ip,iratpt(k,j),iratpt(kp1,j),x,y,z,xh,yh,zh,
 320      . rn(1,j),shifthm)
 321             rcnet(ip) = rcnet(ip)+ signr(j)*5.4548*shifthm
 322    80     continue
 323   100   continue
 324 c
 325 120     continue !PROTON LOOP
 326   116 CONTINUE
 327 C *************BEGIN ANISOTROPY CALCULATION**********************
 328 C ****Locate all backbone C=O and calculate anisotropic shift****
 329         type *,' Calculating CO anisotropy...'
 330         do 30 I=1,iheavyatom
 331           if (aname(I).ne.' CA ') goto 30
 332           DO 35 J=1,15
 333             IF ((ANAME(I+J).EQ.' C  ').AND.(IRES(I).EQ.IRES(I+J)))
 334      1       then
 335              Inext=I+J
 336              GOTO 38
 337             ENDIF
 338 35        CONTINUE
 339           print 95,IRES(I)
 340 95        format (' C atom not found for residue ',I5)
 341         goto 30
 342 38          IF (ANAME(Inext+1).NE.' O  ') then
 343              print 94,IRES(I)
 344 94           format (' O atom not found for residue ',I5)
 345              goto 30
 346             ENDIF
 347 C Atoms CA, C and O now located at I, Inext and (Inext+1)
 348           at3=I
 349           at1=Inext
 350           at2=Inext+1
 351         CALL VEC(at1,at2,at3,vx,vy,vz,r)
 352 C Returns vectors vx,vy,vz, and distance r between at1 and at2
 353         CALL CENTRE(vz,x(at1),y(at1),z(at1),1.1,cen)
 354 C
 355 C Now loop through all H, calculating shift due to this C=O
 356         do 40 K=1,iproton
 357           IF(IRESH(K).EQ.IRES(AT1)) GOTO 40     !SELF SHIFT
 358           CALL RHAT(xh(k),yh(k),zh(k),cen,rh,rlen)
 359 C Returns vector rh, length rlen, between HA and CEN
 360           IF(RLEN.GT.12.0) GOTO 40      !distance cutoff
 361           CALL VPROD(rh,vy,pregam,stheta,tempab)
 362 C Returns sine of angle theta between rh and vy
 363           CALL VPROD(pregam,vx,yy,sgamma,tempab)
 364           calc1=XCO1*((3.0*stheta*stheta)-2.0)
 365           calc2=XCO2*(1.0-(3.0*stheta*stheta*sgamma*sgamma))
 366           calc3=(calc1+calc2)/(3.0*rlen*rlen*rlen)
 367           shift(k)=shift(k)-calc3
 368           anisco(k)=anisco(k)-calc3
 369 40      continue
 370 30      continue
 371 C       if(XCO1.ne.9999.9) goto 9999 !To skip sidechain CO calculation
 372 C *********************************************************
 373 C *******************Sidechain CO from Asn, Gln************
 374         do 530 I=1,iheavyatom
 375         if(res(I).ne.'ASN'.and.res(I).ne.'GLN') goto 530
 376           if (aname(I).ne.' CB ') goto 530
 377           if(res(I).eq.'ASN') then
 378            if (aname(I+1).ne.' CG '.or.aname(I+2).ne.' OD1') then
 379                 print 595, IRES(I)
 380                 goto 530
 381 595       format (' Missing atoms for ASN ',I5)
 382            endif
 383           at1=I+1
 384           at2=I+2
 385           at3=I
 386           else if(res(I).eq.'GLN') then
 387            if (aname(I+1).ne.' CG '.or.aname(I+2).ne.' CD '.or.
 388      1      aname(I+3).ne.' OE1')then
 389                 print 596, IRES(I)
 390                 goto 530
 391 596       format (' Missing atoms for GLN ',I5)
 392            endif
 393           at1=I+2
 394           at2=I+3
 395           at3=I+1
 396          endif
 397         CALL VEC(at1,at2,at3,vx,vy,vz,r)
 398         CALL CENTRE(vz,x(at1),y(at1),z(at1),1.1,cen)
 399 C Now loop through all HA, calculating shift due to this C=O
 400         do 540 K=1,iproton
 401           CALL RHAT(xh(k),yh(k),zh(k),cen,rh,rlen)
 402           IF(RLEN.GT.12.0) GOTO 540
 403           CALL VPROD(rh,vy,pregam,stheta,tempab)
 404           CALL VPROD(pregam,vx,yy,sgamma,tempab)
 405           calc1=XCO1*((3.0*stheta*stheta)-2.0)
 406           calc2=XCO2*(1.0-(3.0*stheta*stheta*sgamma*sgamma))
 407           calc3=(calc1+calc2)/(3.0*rlen*rlen*rlen)
 408           shift(k)=shift(k)-calc3
 409           anisco(k)=anisco(k)-calc3
 410 540     continue
 411 530     continue
 412 C *****************************************************
 413 C *****************Sidechain CO from Asp, Glu**********
 414         do 630 I=1,iheavyatom
 415         if(res(I).ne.'ASP'.and.res(I).ne.'GLU') goto 630
 416           if (aname(I).ne.' CB ') goto 630
 417           if(res(I).eq.'ASP') then
 418            if (aname(I+1).ne.' CG '.or.aname(I+2).ne.' OD1'.or.
 419      1       aname(I+3).ne.' OD2') then
 420                 print 695, IRES(I)
 421                 goto 630
 422 695       format (' Missing atoms for ASP ',I5)
 423            endif
 424           at1=I+1
 425           at2=I+2
 426           at3=I
 427           else if(res(I).eq.'GLU') then
 428            if (aname(I+1).ne.' CG '.or.aname(I+2).ne.' CD '.or.
 429      1      aname(I+3).ne.' OE1'.or.aname(I+4).ne.' OE2')then
 430                 print 696, IRES(I)
 431                 goto 630
 432 696       format (' Missing atoms for GLU ',I5)
 433            endif
 434           at1=I+2
 435           at2=I+3
 436           at3=I+1
 437          endif
 438         do 650 Itmp=1,iproton
 439         calctemp(Itmp)=0.0
 440 650     continue
 441         do 667 Icalc=0,1                !loop for two C-O
 442         at2=at2+Icalc
 443         CALL VEC(at1,at2,at3,vx,vy,vz,r)
 444         CALL CENTRE(vz,x(at1),y(at1),z(at1),1.1,cen)
 445         do 640 K=1,iproton
 446           CALL RHAT(xh(k),yh(k),zh(k),cen,rh,rlen)
 447           IF(RLEN.GT.12.0) GOTO 640
 448           CALL VPROD(rh,vy,pregam,stheta,tempab)
 449           CALL VPROD(pregam,vx,yy,sgamma,tempab)
 450           calc1=XCO1*((3.0*stheta*stheta)-2.0)
 451           calc2=XCO2*(1.0-(3.0*stheta*stheta*sgamma*sgamma))
 452           calc3=(calc1+calc2)/(3.0*rlen*rlen*rlen)
 453         calctemp(K)=calctemp(K)+calc3
 454 640     continue
 455 667     continue
 456           do 664 kk=1,iproton
 457           shift(kk)=shift(kk)-(calctemp(kk)/2.0)
 458           anisco(kk)=anisco(kk)-(calctemp(kk)/2.0)
 459 664       continue
 460 630     continue
 461 9999    continue
 462 C *******NOW DO (O=)C-N TOO *****************************
 463         type *,' Calculating CN anisotropy...'
 464         do 130 I=1,iheavyatom
 465           if (aname(I).ne.' C  ') goto 130
 466             IF (ANAME(I+1).NE.' O  ') then
 467              print 94,IRES(I)
 468              goto 130
 469             ENDIF
 470           do 131 j=1,23
 471             if((aname(I+j).eq.' N  ').and.(ires(I).eq.(ires(i+j)-1)))
 472      1       then
 473              Inext=I+J
 474              GOTO 132
 475             ENDIF
 476 131       CONTINUE
 477              if(ires(I).ne.ires(iheavyatom)) then
 478              print 194,IRES(I)
 479 194          format(' N atom not found after residue ',I5)
 480              endif
 481              goto 130
 482 132       at1=I
 483           at2=Inext
 484           at3=I+1
 485         CALL VEC(at1,at2,at3,vx,vy,vz,r)
 486         rbond=r*0.85    !i.e. 85% along C-N bond
 487         CALL CENTRE(vz,x(at1),y(at1),z(at1),rbond,cen)
 488         do 140 K=1,iproton
 489           CALL RHAT(xh(k),yh(k),zh(k),cen,rh,rlen)
 490           IF(RLEN.GT.12.0) GOTO 140
 491         if(anameh(K).eq.' H  '.and.(iresh(k).eq.ires(Inext)))
 492      $    goto 140
 493           CALL VPROD(rh,vy,pregam,stheta,tempab)
 494           CALL VPROD(pregam,vx,yy,sgamma,tempab)
 495           calc1=XCN1*((3.0*stheta*stheta)-2.0)
 496           calc2=XCN2*(1.0-(3.0*stheta*stheta*sgamma*sgamma))
 497           calc3=(calc1+calc2)/(3.0*rlen*rlen*rlen)
 498           shift(k)=shift(k)-calc3
 499           aniscn(k)=aniscn(k)-calc3
 500 140     continue
 501 130     continue
 502 C *************************************************
 503 C *******Calculate electric field component********
 504 C First find C(alpha)-H(alpha) pair, then work through
 505 C  all C, O and N finding field from them
 506         type *,' Calculating electric field shift...'
 507         do 2000 ie=1,iproton
 508 C skip for NH proton
 509         if(anameh(ie).eq.' H  ')goto 2000
 510         iha=ie
 511         Ez=0.0
 512         do 2010 je=1,iheavyatom
 513 C Find attached heavy atom
 514           if((ires(je).eq.iresh(ie)).and.(aname(je)(3:3).eq.
 515      1      anameh(ie)(3:3))) then
 516           ica=je
 517           goto 2020
 518           endif
 519 2010    continue
 520 2020    continue
 521         call VEC2(iha,ica,vca,rca)
 522 C Returns vector vca and length rca between H(ie) and CA(je)
 523 C Now go through each C and add shift due to this
 524         do 2030 je=1,iheavyatom
 525         if(ires(je).eq.iresh(iha)) goto 2030
 526         if(aname(je).eq.' C  ') then
 527           call VEC2(iha,je,vc,rcc)
 528            if(rcc.gt.6.0) goto 2030     !Ignore for distance>6A
 529           call VSCAL(vca,vc,sc,cthet)
 530           Efact=Qc/(rcc*rcc)
 531           Ez=Ez+(cthet*Efact)
 532         ELSE if(aname(je).eq.' N  ') then
 533           call VEC2(iha,je,vn,rcn)
 534            if(rcn.gt.6.0) goto 2030
 535           call VSCAL(vca,vn,sc,cthet)
 536           Efact=Qn/(rcn*rcn)
 537           Ez=Ez+(cthet*Efact)
 538         ELSE if(aname(je).eq.' O  ') then
 539           call VEC2(iha,je,vo,rco)
 540            if(rco.gt.6.0) goto 2030
 541           call VSCAL(vca,vo,sc,cthet)
 542           Efact=Qo/(rco*rco)
 543           Ez=Ez+(cthet*Efact)
 544         ELSE if(aname(je).eq.' H  ') then
 545           call VEC2(iha,je,vchn,rchn)
 546            if(rchn.gt.6.0) goto 2030
 547           call VSCAL(vca,vchn,sc,cthet)
 548           Efact=Qhn/(rchn*rchn)
 549           Ez=Ez+(cthet*Efact)
 550         endif
 551 2030    continue
 552         shift(ie)=shift(ie)+(Ez*sigmaE)
 553         shiftE(ie)=shiftE(ie)+(Ez*sigmaE)
 554 2000    continue
 555 C ****************END OF CALCULATIONS PROPER********************
 556 C *Correction of Gly HA shift by 0.22 ppm for random coil effect,
 557 C subtract 0.20 from HN shift, subtract 0.65 from HA shift,
 558 C convert to single precision, add random coil shift
 559         do 2040 kkk=1,iproton
 560         shift(kkk)=shift(kkk) + rcnet(kkk)
 561         if(anameh(kkk)(2:3).eq.'HA') shift(kkk)=shift(kkk)-0.65
 562         if(anameh(kkk)(2:3).eq.'H ') shift(kkk)=shift(kkk)-0.20
 563         sum(kkk)=rcnet(kkk)+anisco(kkk)+aniscn(kkk)+shiftE(kkk)
 564         if((resh(kkk).eq.'GLY').and.(anameh(kkk).eq.'1HA '.or.
 565      1   anameh(kkk).eq.'2HA '.or.anameh(kkk).eq.' HA1'.or.
 566      2   anameh(kkk).eq.' HA2')) shift(kkk)=shift(kkk)+0.22
 567           do 2050 L=1,179
 568            IF ((datres(L).eq.resh(kkk)).and.(datnam(L)(2:4).eq.
 569      1     anameh(kkk)(2:4))) THEN
 570 C This next bit replaces HA shifts for aromatics by 4.45 ppm
 571 C  (gives roughly best fit to data)
 572             if((resh(kkk).eq.'TYR'.or.resh(kkk).eq.'PHE'.or.
 573      1      resh(kkk).eq.'TRP'.or.resh(kkk).eq.'HIS').and.anameh(kkk).
 574      2      eq.' HA ') then
 575               final(kkk)=shift(kkk)+4.45
 576               shift(kkk)=shift(kkk)+4.45-datshift(L)
 577             ELSE
 578               final(kkk)=shift(kkk) + datshift(L)
 579             ENDIF
 580           ENDIF
 581 2050      continue
 582 2040    continue
 583 9990    continue
 584
 585 C ***************** Output SHIFT to channel 3**********************
 586         if(yn.ne.'Y'.and.yn.ne.'y') goto 2100
 587         do 909 ii=1,iprot
 588         DO 910 KK=1,iproton
 589         IF(IEXP(ii).EQ.IRESH(KK).AND.pexp(ii).eq.prott.and.
 590      1   atexp(ii)(1:2).eq.anameh(kk)(2:3)) THEN
 591           exptln(kk)=eexp(ii)
 592 C exptln is the experimental shift: now subtract random coil
 593           do 2051 L=1,179
 594            IF ((datres(L).eq.resh(kk)).and.(datnam(L)(1:3).eq.
 595      1     anameh(kk)(1:3))) then
 596            exptln(kk)=exptln(kk) - datshift(L)
 597            ENDIF
 598 2051      continue
 599           exptl(kk)=exptln(kk)
 600           goto 909
 601         ENDIF
 602 910     continue
 603 909     continue
 604         if(imodel.eq.1) goto 501
 605         write (3,988)
 606 988     format(' Proton name, followed by calc. and observed shift and
 607      . difference')
 608         write (3,2052)
 609 2052    format(' (as shift - random coil shift)')
 610         do 500 iprnt=1,iproton
 611           if(exptln(iprnt).ne.0.00)
 612      1     write(3,991) iprnt,iresh(iprnt),resh(iprnt),anameh(iprnt),
 613      1     shift(iprnt),exptl(iprnt),(shift(iprnt)-exptl(iprnt))
 614 991     format(I5,I5,1X,A3,2X,A4,3F10.6)
 615 500     continue
 616 501     continue
 617         if(imodel.eq.1) then
 618         write(3,775) nmodel
 619 775     format(' Result for model number',I8)
 620         endif
 621         CALL STATS(shift,EXPTL,iproton,anameh)
 622         type *,'Do you want output sorted? [N]'
 623         read(5,970)ans
 624         if(ans.ne.'Y'.and.ans.ne.'y') goto 1000
 625         itemp=0
 626         do 505 I=1,iproton
 627         if(exptl(I).eq.0.0) goto 505
 628         itemp=itemp+1
 629         shiftt(itemp)=shift(I)
 630         exptlt(itemp)=exptl(I)
 631 505     continue
 632         CALL INDEXX(itemp,shiftt,indx)
 633         write (3,990)
 634 C990    format('    Calc       Exptl       Diff')
 635 990     format('    Calc       Diff')
 636         do 506 I=1,itemp
 637 C       write (3,989) shiftt(indx(I)),exptlt(indx(I)),
 638         write (3,983) shiftt(indx(I)),
 639      1   (shiftt(indx(I))-exptlt(indx(I)))
 640 989     format(3F10.4)
 641 983     format(2F10.4)
 642 506     continue
 643         goto 1000
 644 C **************Normal output starts here***************************
 645 2100    continue
 646 C       write (3,987) FN1
 647 987     format(' Shift calculated for protons from ',A)
 648 C       write (3,986)
 649 986     format(' Added 0.22 for Gly, subtracted 0.65 (HA only), added
 650      1random coil shift')
 651 C       write (3,985)
 652 985     format('           Proton       ring     anisCO    anisCN
 653      1sigmaE   sum       final')
 654         do 510 ip=1,iproton
 655         if(icalculate.eq.1.and.anameh(ip)(2:3).ne.'HA') goto 510
 656         if(icalculate.eq.2.and.anameh(ip)(2:3).ne.'H ') goto 510
 657 C Now do averaging for Phe,Tyr,methyls
 658         if((resh(ip).eq.'VAL'.and.(anameh(ip).eq.'1HG1'.or.anameh(ip).
 659      1eq.'1HG2')).or.(resh(ip).eq.'LEU'.and.(anameh(ip).eq.'1HD1'.or.
 660      2anameh(ip).eq.'1HD2')).or.(resh(ip).eq.'ILE'.and.(anameh(ip).eq.
 661      3'1HG2'.or.anameh(ip).eq.'1HD1')).or.(resh(ip).eq.'ALA'.and.
 662      4anameh(ip).eq.'1HB ').or.(resh(ip).eq.'THR'.and.anameh(ip).eq.
 663      5'1HG2')) THEN
 664           final(ip)=final(ip)-sum(ip)
 665           sum(ip)=(sum(ip)+sum(ip+1)+sum(ip+2))/3.0
 666           anisco(ip)=(anisco(ip)+anisco(ip+1)+anisco(ip+2))/3.0
 667           aniscn(ip)=(aniscn(ip)+aniscn(ip+1)+aniscn(ip+2))/3.0
 668           shiftE(ip)=(shiftE(ip)+shiftE(ip+1)+shiftE(ip+2))/3.0
 669           shift(ip)=(shift(ip)+shift(ip+1)+shift(ip+2))/3.0
 670           final(ip)=final(ip)+sum(ip)
 671           final(ip+1)=0.0
 672           final(ip+2)=0.0
 673         endif
 674         if(resh(ip).eq.'TYR'.and.anameh(ip).eq.' HD1') then
 675           if(anameh(ip+1).eq.' HD2') then
 676            final(ip)=final(ip)-sum(ip)
 677            sum(ip)=(sum(ip)+sum(ip+1))/2.0
 678            anisco(ip)=(anisco(ip)+anisco(ip+1))/2.0
 679            aniscn(ip)=(aniscn(ip)+aniscn(ip+1))/2.0
 680            shiftE(ip)=(shiftE(ip)+shiftE(ip+1))/2.0
 681            shift(ip)=(shift(ip)+shift(ip+1))/2.0
 682            final(ip)=final(ip)+sum(ip)
 683            final(ip+1)=0.0
 684           else
 685            type 940,resh(ip),iresh(ip)
 686 940       format(' Ring protons in unexpected order for ',A3,I4)
 687           endif
 688         endif
 689         if(resh(ip).eq.'TYR'.and.anameh(ip).eq.' HE1') then
 690           if(anameh(ip+1).eq.' HE2') then
 691            final(ip)=final(ip)-sum(ip)
 692            sum(ip)=(sum(ip)+sum(ip+1))/2.0
 693            anisco(ip)=(anisco(ip)+anisco(ip+1))/2.0
 694            aniscn(ip)=(aniscn(ip)+aniscn(ip+1))/2.0
 695            shiftE(ip)=(shiftE(ip)+shiftE(ip+1))/2.0
 696            shift(ip)=(shift(ip)+shift(ip+1))/2.0
 697            final(ip)=final(ip)+sum(ip)
 698            final(ip+1)=0.0
 699           else
 700            type 940,resh(ip),iresh(ip)
 701           endif
 702         endif
 703         if(resh(ip).eq.'PHE'.and.anameh(ip).eq.' HD1') then
 704           if(anameh(ip+1).eq.' HD2') then
 705            final(ip)=final(ip)-sum(ip)
 706            sum(ip)=(sum(ip)+sum(ip+1))/2.0
 707            anisco(ip)=(anisco(ip)+anisco(ip+1))/2.0
 708            aniscn(ip)=(aniscn(ip)+aniscn(ip+1))/2.0
 709            shiftE(ip)=(shiftE(ip)+shiftE(ip+1))/2.0
 710            shift(ip)=(shift(ip)+shift(ip+1))/2.0
 711            final(ip)=final(ip)+sum(ip)
 712            final(ip+1)=0.0
 713           else
 714            type 940,resh(ip),iresh(ip)
 715           endif
 716         endif
 717         if(resh(ip).eq.'PHE'.and.anameh(ip).eq.' HE1') then
 718           if(anameh(ip+1).eq.' HE2') then
 719            final(ip)=final(ip)-sum(ip)
 720            sum(ip)=(sum(ip)+sum(ip+1))/2.0
 721            anisco(ip)=(anisco(ip)+anisco(ip+1))/2.0
 722            aniscn(ip)=(aniscn(ip)+aniscn(ip+1))/2.0
 723            shiftE(ip)=(shiftE(ip)+shiftE(ip+1))/2.0
 724            shift(ip)=(shift(ip)+shift(ip+1))/2.0
 725            final(ip)=final(ip)+sum(ip)
 726            final(ip+1)=0.0
 727           else
 728            type 940,resh(ip),iresh(ip)
 729           endif
 730         endif
 731 C ******Write out results*********************
 732         if(final(ip).ne.0.0) write(3,984)ip,iresh(ip),resh(ip),
 733      1   anameh(ip),rcnet(ip),anisco(ip),aniscn(ip),shiftE(ip),
 734      2   shift(ip),final(ip)
 735 984     format(I5,I5,1X,A3,2X,A4,6F10.5)
 736 510     CONTINUE
 737 1000    continue
 738         imodel=1 !to make it work for single structure
 739         goto 776
 740 C       Hacked DvdS 10/98
 741  777    continue
 742         print *,'Closing unit 3'
 743         close(3)
 744 C
 745         STOP
 746         END
 747 C *******************************************************
 748         SUBROUTINE VEC(at1,at2,at3,vx,vy,vz,r)
 749         COMMON X,Y,Z
 750         DIMENSION X(5000),Y(5000),Z(5000),vx(3),vy(3),vz(3)
 751         DIMENSION vvx(3),vvy(3),vvz(3)
 752         INTEGER at1,at2,at3
 753         x1=x(at1)
 754         y1=y(at1)
 755         z1=z(at1)
 756         x2=x(at2)
 757         y2=y(at2)
 758         z2=z(at2)
 759         x3=x(at3)
 760         y3=y(at3)
 761         z3=z(at3)
 762         vvz(1)=x2-x1
 763         vvz(2)=y2-y1
 764         vvz(3)=z2-z1
 765         CALL UNITV(vvz,vz,r)
 766         vvz(1)=x3-x1
 767         vvz(2)=y3-y1
 768         vvz(3)=z3-z1
 769         CALL VPROD(vz,vvz,vvy,sthet,tempab)
 770         vvy(1)=tempab
 771         CALL UNITV(vvy,vy,D)
 772         CALL VPROD(vz,vy,vvx,sthet,tempab)
 773         CALL UNITV(vvx,vx,D)
 774         return
 775         end
 776 C ******************************************************
 777         SUBROUTINE UNITV(U,U1,D)
 778         DIMENSION U(3),U1(3)
 779         D=0.0
 780         DO 1 JJ=1,3
 781         D=D+(U(JJ)*U(JJ))
 782 1       CONTINUE
 783         D=SQRT(D)
 784         If (D.EQ.0.0) type *,' D is zero in UNITV'
 785         DO 2 JJ=1,3
 786         U1(JJ)=U(JJ)/D
 787 2       CONTINUE
 788         RETURN
 789         END
 790 C **********************************************************
 791         SUBROUTINE VPROD(A,B,AB,STHET,tempab)
 792         DIMENSION A(3),B(3),AB(3)
 793         AB(2)=A(3)*B(1) - A(1)*B(3)
 794         AB(1)=A(2)*B(3) - A(3)*B(2)
 795         AB(3)=A(1)*B(2) - A(2)*B(1)
 796         tempab=AB(1)
 797 c This is a stupid thing to do, but it's the only way I can get it
 798 c  to behave!
 799         R1=SQRT(AB(1)*AB(1) + AB(2)*AB(2) + AB(3)*AB(3))
 800         RA=SQRT(A(1)*A(1) + A(2)*A(2) + A(3)*A(3))
 801         RB=SQRT(B(1)*B(1) + B(2)*B(2) + B(3)*B(3))
 802           IF(RA.EQ.0.0.OR.RB.EQ.0.0) THEN
 803             type *,' VPROD Zero divide...ignore',RA,RB
 804             STHET=0.0
 805           ELSE
 806             STHET=R1/(RA*RB)
 807           ENDIF
 808         RETURN
 809         END
 810 C *******************************************************
 811         SUBROUTINE CENTRE(vz,a1,a2,a3,dist,cen)
 812         DIMENSION vz(3),cen(3)
 813         cen(1)=a1+(vz(1)*dist)
 814         cen(2)=a2+(vz(2)*dist)
 815         cen(3)=a3+(vz(3)*dist)
 816         return
 817         end
 818 C ********************************************************
 819         SUBROUTINE RHAT(ax,ay,az,cen,rh,rlen)
 820         DIMENSION cen(3),rh(3)
 821         rh(1)=ax-cen(1)
 822         rh(2)=ay-cen(2)
 823         rh(3)=az-cen(3)
 824         rlen=sqrt(rh(1)*rh(1)+rh(2)*rh(2)+rh(3)*rh(3))
 825         return
 826         end
 827 C ********************************************************
 828         SUBROUTINE STATS(SSHIFT,EXPTL,iproton,anameh)
 829         dimension sshift(3000),exptl(3000),anameh(3000)
 830         dimension zshift(3000),zexptl(3000)
 831 C Calculates mean and sd of iproton values of shift and exptl
 832 C  also mean and sd of (shift-exptl) and regression coeff
 833 C Values of exptl of 0.00 presumably not found in protha.out
 834 C and are excluded
 835         itemp=0
 836         do 5 I=1,iproton
 837         if(exptl(I).eq.0.0) goto 5
 838         itemp=itemp+1
 839         if(anameh(I).eq.'1HA '.and.anameh(I+1).eq.'2HA ') then
 840           sshift(I)=(sshift(I)+sshift(I+1))/2.0
 841         endif
 842         zshift(itemp)=sshift(I)
 843         zexptl(itemp)=exptl(I)
 844 5       continue
 845         sm1=0.0
 846         sm2=0.0
 847         sms1=0.0
 848         sms2=0.0
 849         sm3=0.0
 850         sms3=0.0
 851         sm4=0.0
 852 C sm1,sm2 are accumulated sums of shift and exptl
 853 C sms1,sms2 are accumulated sums of shift**2 and exptl**2
 854         do 600 L=1,itemp
 855         sm1=sm1+zshift(L)
 856         sm2=sm2+zexptl(L)
 857         sm3=sm3+(zshift(L)-zexptl(L))
 858 600     continue
 859         sm1=sm1/float(itemp)
 860         sm2=sm2/float(itemp)
 861         sm3=sm3/float(itemp)
 862         do 605 L=1,itemp
 863         sms1=sms1+((zshift(L)-sm1)*(zshift(L)-sm1))
 864         sms2=sms2+((zexptl(L)-sm2)*(zexptl(L)-sm2))
 865         sms3=sms3+((zshift(L)-zexptl(L)-sm3)*(zshift(L)-zexptl(L)-sm3))
 866 605     continue
 867         sms1=sqrt(sms1/float(itemp))
 868         sms2=sqrt(sms2/float(itemp))
 869         sms3=sqrt(sms3/float(itemp))
 870         do 610 L=1,itemp
 871         sm4=sm4+((zshift(L)-sm1)*(zexptl(L)-sm2))
 872 610     continue
 873         sm4=sm4/(float(itemp)*sms1*sms2)
 874         write (3,1000) sm1,sm2
 875 1000    format(' Means of calc and exptl shifts: ',F6.4,
 876      1   1X,F6.4)
 877         write (3,1001) sms1,sms2
 878 1001    format(' SD of calc and exptl shifts: ',F6.4,
 879      1   1X,F6.4)
 880         write (3,1002) sm3,sms3
 881 1002    format(' Mean and SD of (calc-exptl shift): ',F6.4,
 882      1   1X,F6.4)
 883         write (3,1003) sm4,itemp
 884 1003    format (' Correlation coefficient calc vs exptl: ',F7.5,
 885      1   ' (',I4,' protons )')
 886         return
 887         end
 888 C ************************************************
 889         SUBROUTINE VEC2(iha,ica,vca,rca)
 890         DIMENSION vca(3),x(5000),y(5000),z(5000)
 891         DIMENSION xh(3000),yh(3000),zh(3000)
 892         COMMON x,y,z,xh,yh,zh
 893         vca(1)=xh(iha)-x(ica)
 894         vca(2)=yh(iha)-y(ica)
 895         vca(3)=zh(iha)-z(ica)
 896         rca=sqrt((vca(1)*vca(1))+(vca(2)*vca(2))+(vca(3)*vca(3)))
 897         return
 898         end
 899 C *******************************************
 900        SUBROUTINE VSCAL(A,B,SC,CTHET)
 901        DIMENSION A(3),B(3)
 902        SC=A(1)*B(1)+A(2)*B(2)+A(3)*B(3)
 903        RA=SQRT(A(1)*A(1)+A(2)*A(2)+A(3)*A(3))
 904        RB=SQRT(B(1)*B(1)+B(2)*B(2)+B(3)*B(3))
 905        CTHET=SC/(RA*RB)
 906        RETURN
 907        END
 908 C **********************************************
 909         SUBROUTINE INDEXX(N,ARRIN,INDX)
 910 C Indexes an array ARRIN of length N, i.e. outputs array INDX
 911 C such that ARRIN(INDX(J)) is in ascending order.
 912 C Taken from 'Numerical Recipes', W.H.Press et al, CUP 1989
 913         DIMENSION ARRIN(3000),INDX(3000)
 914         DO 11 J=1,N
 915                 INDX(J)=J
 916 11              CONTINUE
 917         IF(N.EQ.1) RETURN
 918         L=N/2+1
 919         IR=N
 920 10      CONTINUE
 921            IF(L.GT.1) THEN
 922                 L=L-1
 923                 INDXT=INDX(L)
 924                 Q=ARRIN(INDXT)
 925            ELSE
 926                 INDXT=INDX(IR)
 927                 Q=ARRIN(INDXT)
 928                 INDX(IR)=INDX(1)
 929                 IR=IR-1
 930                 IF(IR.EQ.1) THEN
 931                   INDX(1)=INDXT
 932                   RETURN
 933                 ENDIF
 934            ENDIF
 935            I=L
 936            J=L+L
 937 20         IF(J.LE.IR) THEN
 938                 IF(J.LT.IR) THEN
 939                   IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
 940                 ENDIF
 941                 IF(Q.LT.ARRIN(INDX(J))) THEN
 942                   INDX(I)=INDX(J)
 943                   I=J
 944                   J=J+J
 945                 ELSE
 946                   J=IR+1
 947                 ENDIF
 948            GO TO 20
 949            ENDIF
 950            INDX(I)=INDXT
 951         GOTO 10
 952         END
 953 c**********************************************************
 954       subroutine hm(ip,ir1,ir2,x,y,z,xh,yh,zh,rn,shhm)
 955 c
 956 c  Subroutine Haigh-Mallion:
 957 c  -- given proton ip, and ring atoms ir1 and ir2, with
 958 c        coordinates x and ring normal vector rn;
 959 c  -- return Haigh-Maillion shift contribution shhm
 960 c
 961         dimension x(5000),y(5000),z(5000),xh(3000),yh(3000)
 962       dimension rn(3),r1(3),r2(3),zh(3000)
 963 c
 964 c --- extract coordinates from array x,y,z
 965 c
 966         r1(1)=x(ir1) - xh(ip)
 967         r1(2)=y(ir1) - yh(ip)
 968         r1(3)=z(ir1) - zh(ip)
 969         r2(1)=x(ir2) - xh(ip)
 970         r2(2)=y(ir2) - yh(ip)
 971         r2(3)=z(ir2) - zh(ip)
 972 c
 973 c  -- compute triple scalar product: later versions could
 974 c       hard-code this for efficiency, but this form is
 975 c       easier to read. (Triple scalar product=vol. of
 976 c       parallelepiped, =area of desired triangle)
 977 c
 978       s12 = r1(1)*(r2(2)*rn(3) - r2(3)*rn(2))
 979      .    + r1(2)*(r2(3)*rn(1) - r2(1)*rn(3))
 980      .    + r1(3)*(r2(1)*rn(2) - r2(2)*rn(1))
 981 c
 982 c  -- get radial factors
 983 c
 984       r1sq = r1(1)*r1(1) + r1(2)*r1(2) + r1(3)*r1(3)
 985       r2sq = r2(1)*r2(1) + r2(2)*r2(2) + r2(3)*r2(3)
 986       if (r1sq.eq.0.0 .or. r2sq.eq.0.0) then
 987         write(6,*) 'Geometry error in hm: ',ip,r1,r2
 988         stop
 989       end if
 990       temp = (1./r1sq**1.5 + 1./r2sq**1.5)*0.5  !distance bit
 991       shhm = s12*temp
 992       return
 993       end
 994 C****************************************************************
 995       subroutine plane (i1, i2, i3, x,y,z, rn, cent)
 996 c**************************************************************
 997 c  --- given three atoms i1,i2 and i3, and coordinates x,y,z
 998 c    - returns rn(i) [i=1,2,3] components of the normalized
 999 c      vector normal to the plane containing the three
1000 c      points and cent(1,2,3), the centre of the three atom.
1001 c
1002       dimension x(5000),y(5000),z(5000)
1003       dimension rn(3),cent(3)
1004 c
1005         x1 = x(i1)
1006         y1 = y(i1)
1007         z1 = z(i1)
1008         x2 = x(i2)
1009         y2 = y(i2)
1010         z2 = z(i2)
1011         x3 = x(i3)
1012         y3 = y(i3)
1013         z3 = z(i3)
1014 c
1015 c       ----- coefficients of the equation for the plane of atoms 1-3
1016 c
1017         ax = y1*z2 - y2*z1 + y3*z1 - y1*z3 + y2*z3 - y3*z2
1018         ay = -(x1*z2 - x2*z1 + x3*z1 - x1*z3 + x2*z3 - x3*z2)
1019         az = x1*y2 - x2*y1 + x3*y1 - x1*y3 + x2*y3 - x3*y2
1020         anorm = 1./sqrt(ax*ax + ay*ay + az*az)
1021 c
1022 c       ----- normalize to standard form for plane equation (i.e. such
1023 c       ----- that length of the vector "a" is unity
1024 c
1025         rn(1) = ax*anorm        !ring normal unit vector
1026         rn(2) = ay*anorm        !=(2-1)x(3-1)
1027         rn(3) = az*anorm        !divided by its magnitude
1028 c
1029         cent(1) = (x1+x2+x3)/3.
1030         cent(2) = (y1+y2+y3)/3.
1031         cent(3) = (z1+z2+z3)/3.
1032 c
1033         return
1034         end
1035 C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1036         subroutine addprot(x,y,z,xh,yh,zh,Iheavyatom,Iproton,
1037      &   aname,anameh,res,resh,ires,iresh)
1038 C Does not add most exchangeable protons: N-terminal NH3,
1039 C  Glu,Gln,Asp,Asn sidechain, Arg sidechain,Tyr HH, His HD1,
1040 C  Trp HE1, Ser & Thr OH, Cys SH, Lys NH3
1041 C Goes through heavy atom list and adds in simple geometric way
1042 C NB Looks in neighbourhood of attached heavy atom for bonding
1043 C atoms. If the atom ordering is odd, this could cause problems.
1044 C Usually this would mean that the proton does not get added;
1045 C very occasionally it could lead to a proton being put on wrong.
1046         DIMENSION x(5000),y(5000),z(5000),xh(3000),yh(3000),zh(3000)
1047         DIMENSION aname(5000),res(5000),ires(5000)
1048         DIMENSION anameh(3000),resh(3000),iresh(3000)
1049         CHARACTER aname*4,anameh*4,res*3,resh*3
1050         iproton=0
1051 C NB This effectively deletes any protons that already were in the file
1052         do 10 I=1,iheavyatom
1053         iat2=0
1054         iat3=0
1055         iat4=0
1056         if(aname(I).eq.' N  '.and.ires(I).gt.1.and.res(I).
1057      &    ne.'PRO') then
1058          iat1=I
1059          do 11 J=1,13
1060          if(aname(I-J).eq.' C  '.and.ires(I-J).eq.ires(I)-1) then
1061          iat2=I-J
1062          goto 12
1063          endif
1064 11       continue
1065 12       continue
1066          do 13 J=1,15
1067          if(aname(I+J).eq.' CA '.and.ires(I+J).eq.ires(I)) then
1068          iat3=I+J
1069          goto 14
1070          endif
1071 13       continue
1072 14       continue
1073          if(iat2.gt.0.and.iat3.gt.0) then
1074          iproton=iproton+1
1075          iheavyatom=iheavyatom+1
1076          call addsp2(x,y,z,xh,yh,zh,iproton,iat1,iat2,iat3,1)
1077          anameh(iproton)=' H  '
1078          resh(iproton)=res(iat1)
1079          iresh(iproton)=ires(iat1)
1080          x(iheavyatom)=xh(iproton)
1081          y(iheavyatom)=yh(iproton)
1082          z(iheavyatom)=zh(iproton)
1083          aname(iheavyatom)=' H  '
1084          res(iheavyatom)=res(iat1)
1085          ires(iheavyatom)=ires(iat1)
1086 C The final 1 signals that it is an NH ie use NH bond distance
1087          endif
1088         else if (aname(I).eq.' CA '.and.res(I).ne.'GLY') then
1089          iat1=I
1090          do 21 J=-5,3
1091          if(aname(I+J).eq.' N  '.and.ires(I+J).eq.ires(I)) then
1092          iat2=I+J
1093          goto 22
1094          endif
1095 21       continue
1096 22       continue
1097          do 23 J=1,15
1098          if(aname(I+J).eq.' C  '.and.ires(I+J).eq.ires(I)) then
1099          iat3=I+J
1100          goto 24
1101          endif
1102 23       continue
1103 24       continue
1104          do 25 J=1,15
1105          if(aname(I+J).eq.' CB '.and.ires(I+J).eq.ires(I)) then
1106          iat4=I+J
1107          goto 26
1108          endif
1109 25       continue
1110 26       continue
1111          if(iat2.gt.0.and.iat3.gt.0.and.iat4.gt.0) then
1112          iproton=iproton+1
1113          call addonesp3(x,y,z,xh,yh,zh,iproton,iat1,iat2,iat3,iat4)
1114          anameh(iproton)=' HA '
1115          resh(iproton)=res(iat1)
1116          iresh(iproton)=ires(iat1)
1117          endif
1118         else if((aname(I).eq.' CB '.and.(res(I).eq.'THR'.or.
1119      &   res(I).eq.'VAL'.or.res(I).eq.'ILE')).or.(aname(I).
1120      &   eq.' CG '.and.res(I).eq.'LEU')) then
1121          iat1=I
1122          do 31 J=1,6
1123          if(aname(I-J).eq.' CA '.and.ires(I-J).eq.ires(I).
1124      &    and.res(I).ne.'LEU') then
1125          iat2=I-J
1126          goto 32
1127          endif
1128 31       continue
1129 32       continue
1130          do 33 J=1,6
1131          if(res(I-J).eq.'LEU'.and.aname(I-J).eq.' CB ') then
1132          iat2=I-J
1133          goto 34
1134          endif
1135 33       continue
1136 34       continue
1137          do 35 J=1,3
1138          if(res(I+J).eq.'THR'.and.aname(I+J)(2:3).eq.'OG')
1139      $    iat3=I+J
1140          if(res(I+J).eq.'THR'.and.aname(I+J)(2:3).eq.'CG')
1141      $    iat4=I+J
1142          if(res(I+J).eq.'VAL'.and.aname(I+J).eq.' CG1') iat3=I+J
1143          if(res(I+J).eq.'VAL'.and.aname(I+J).eq.' CG2') iat4=I+J
1144          if(res(I+J).eq.'ILE'.and.aname(I+J).eq.' CG1') iat3=I+J
1145          if(res(I+J).eq.'ILE'.and.aname(I+J).eq.' CG2') iat4=I+J
1146          if(res(I+J).eq.'LEU'.and.aname(I+J).eq.' CD1') iat3=I+J
1147          if(res(I+J).eq.'LEU'.and.aname(I+J).eq.' CD2') iat4=I+J
1148 35       continue
1149          if(iat2.gt.0.and.iat3.gt.0.and.iat4.gt.0) then
1150          iproton=iproton+1
1151          call addonesp3(x,y,z,xh,yh,zh,iproton,iat1,iat2,iat3,iat4)
1152          if(res(I).eq.'LEU') then
1153           anameh(iproton)=' HG '
1154          else
1155           anameh(iproton)=' HB '
1156          endif
1157          resh(iproton)=res(iat1)
1158          iresh(iproton)=ires(iat1)
1159          endif
1160         else if(aname(I).eq.' CA '.and.res(I).eq.'GLY') then
1161          iat1=I
1162          do 41 J=-2,2
1163          if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' N  ')iat3=I+J
1164          if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' C  ')iat2=I+J
1165 41       continue
1166          if(iat2.gt.0.and.iat3.gt.0) then
1167          iproton=iproton+2
1168          call addtwosp3(x,y,z,xh,yh,zh,iproton,iat1,iat2,iat3)
1169          anameh(iproton-1)='1HA '
1170          anameh(iproton)='2HA '
1171          resh(iproton-1)=res(iat1)
1172          resh(iproton)=res(iat1)
1173          iresh(iproton-1)=ires(iat1)
1174          iresh(iproton)=ires(iat1)
1175          endif
1176         else if(aname(I).eq.' CB '.and.(res(I).eq.'CYS'.or.res(I).
1177      1   eq.'ASP'.or.res(I).eq.'GLU'.or.res(I).eq.'PHE'.or.res(I).
1178      2   eq.'HIS'.or.res(I).eq.'LYS'.or.res(I).eq.'LEU'.or.res(I).
1179      3   eq.'MET'.or.res(I).eq.'ASN'.or.res(I).eq.'PRO'.or.res(I).
1180      4   eq.'GLN'.or.res(I).eq.'ARG'.or.res(I).eq.'SER'.or.res(I).
1181      5   eq.'TRP'.or.res(I).eq.'TYR')) then
1182          iat1=I
1183          do 42 J=-3,3
1184          if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CA ')iat2=I+J
1185          if(ires(I+J).eq.ires(I).and.aname(I+J)(3:3).eq.'G')
1186      &    iat3=I+J
1187 42      continue
1188          if(iat2.gt.0.and.iat3.gt.0) then
1189          iproton=iproton+2
1190          call addtwosp3(x,y,z,xh,yh,zh,iproton,iat1,iat2,iat3)
1191          anameh(iproton-1)='1HB '
1192          anameh(iproton)='2HB '
1193          resh(iproton-1)=res(iat1)
1194          resh(iproton)=res(iat1)
1195          iresh(iproton-1)=ires(iat1)
1196          iresh(iproton)=ires(iat1)
1197          endif
1198         else if((aname(I).eq.' CG '.and.(res(I).eq.'GLU'.or.res(I).
1199      1   eq.'LYS'.or.res(I).eq.'MET'.or.res(I).eq.'PRO'.or.res(I).
1200      2   eq.'GLN'.or.res(I).eq.'ARG')).or.(aname(I).eq.' CG1'.and.
1201      3   res(I).eq.'ILE')) then
1202          iat1=I
1203          do 43 J=-3,3
1204          if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CB ')iat2=I+J
1205          if(ires(I+J).eq.ires(I).and.aname(I+J)(3:3).eq.'D')
1206      &    iat3=I+J
1207 43      continue
1208          if(iat2.gt.0.and.iat3.gt.0) then
1209          iproton=iproton+2
1210          call addtwosp3(x,y,z,xh,yh,zh,iproton,iat1,iat2,iat3)
1211          anameh(iproton-1)='1HG '
1212          anameh(iproton)='2HG '
1213          if(res(I).eq.'ILE')anameh(iproton-1)='1HG1'
1214          if(res(I).eq.'ILE')anameh(iproton)='2HG1'
1215          resh(iproton-1)=res(iat1)
1216          resh(iproton)=res(iat1)
1217          iresh(iproton-1)=ires(iat1)
1218          iresh(iproton)=ires(iat1)
1219          endif
1220         else if(aname(I).eq.' CD '.and.(res(I).eq.'LYS'.or.res(I).
1221      1   eq.'ARG'.or.res(I).eq.'PRO')) then
1222          iat1=I
1223          do 44 J=-6,3
1224          if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CG ')iat2=I+J
1225          if(ires(I+J).eq.ires(I).and.aname(I+J)(3:3).eq.'E')
1226      &   iat3=I+J
1227          if(ires(I+J).eq.ires(I).and.res(I).eq.'PRO'.and.
1228      &    aname(I+J).eq.' N  ')iat3=I+J
1229 44       continue
1230          if(iat2.gt.0.and.iat3.gt.0) then
1231          iproton=iproton+2
1232          call addtwosp3(x,y,z,xh,yh,zh,iproton,iat1,iat2,iat3)
1233          anameh(iproton-1)='1HD '
1234          anameh(iproton)='2HD '
1235          resh(iproton-1)=res(iat1)
1236          resh(iproton)=res(iat1)
1237          iresh(iproton-1)=ires(iat1)
1238          iresh(iproton)=ires(iat1)
1239          endif
1240         else if(aname(I).eq.' CE '.and.res(I).eq.'LYS') then
1241          iat1=I
1242          do 45 J=-3,3
1243          if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CD ')iat2=I+J
1244          if(ires(I+J).eq.ires(I).and.aname(I+J)(3:3).eq.'Z')
1245      &   iat3=I+J
1246 45       continue
1247          if(iat2.gt.0.and.iat3.gt.0) then
1248          iproton=iproton+2
1249          call addtwosp3(x,y,z,xh,yh,zh,iproton,iat1,iat2,iat3)
1250          anameh(iproton-1)='1HE '
1251          anameh(iproton)='2HE '
1252          resh(iproton-1)=res(iat1)
1253          resh(iproton)=res(iat1)
1254          iresh(iproton-1)=ires(iat1)
1255          iresh(iproton)=ires(iat1)
1256          endif
1257         else if(aname(I).eq.' CB '.and.res(I).eq.'ALA') then
1258          iat1=I
1259          do 51 J=1,5
1260          if(ires(I-J).eq.ires(I).and.aname(I-J).eq.' CA ')iat2=I-J
1261          if(ires(I-J).eq.ires(I).and.aname(I-J).eq.' N  ')iat3=I-J
1262 51       continue
1263          if(iat2.gt.0.and.iat3.gt.0)goto 52
1264        else if((aname(I)(2:3).eq.'CG'.and.(res(I).eq.'VAL'.or.res(I).
1265      1  eq.'THR')).or.(aname(I).eq.' CG2'.and.res(I).eq.'ILE'))then
1266         iat1=I
1267         do 53 J=1,5
1268         if(ires(I-J).eq.ires(I).and.aname(I-J).eq.' CB ')iat2=I-J
1269         if(ires(I-J).eq.ires(I).and.aname(I-J).eq.' CA ')iat3=I-J
1270 53      continue
1271         if(iat2.gt.0.and.iat3.gt.0)goto 52
1272        else if(aname(I)(2:3).eq.'CD'.and.(res(I).eq.'LEU'.or.res(I).
1273      1   eq.'ILE'))then
1274         iat1=I
1275         do 54 J=1,5
1276         if(ires(I-J).eq.ires(I).and.aname(I-J)(2:3).eq.'CG')iat2=I-J
1277         if(ires(I-J).eq.ires(I).and.aname(I-J).eq.' CB ')iat3=I-J
1278 54      continue
1279         if(iat2.gt.0.and.iat3.gt.0)goto 52
1280        else if(aname(I)(2:3).eq.'CE'.and.res(I).eq.'MET')then
1281         iat1=I
1282         do 55 J=1,5
1283         if(ires(I-J).eq.ires(I).and.aname(I-J).eq.' SD ')iat2=I-J
1284         if(ires(I-J).eq.ires(I).and.aname(I-J).eq.' CG ')iat3=I-J
1285 55      continue
1286         if(iat2.gt.0.and.iat3.gt.0) then
1287 52      iproton=iproton+3
1288         call addthreesp3(x,y,z,xh,yh,zh,iproton,iat1,iat2,iat3)
1289         anameh(iproton-2)='1H  '
1290         anameh(iproton-1)='2H  '
1291         anameh(iproton)='3H  '
1292         anameh(iproton-2)(3:4)=aname(iat1)(3:4)
1293         anameh(iproton-1)(3:4)=aname(iat1)(3:4)
1294         anameh(iproton)(3:4)=aname(iat1)(3:4)
1295         resh(iproton-2)=res(iat1)
1296         resh(iproton-1)=res(iat1)
1297         resh(iproton)=res(iat1)
1298         iresh(iproton-2)=ires(iat1)
1299         iresh(iproton-1)=ires(iat1)
1300         iresh(iproton)=ires(iat1)
1301         endif
1302        else if((res(I).eq.'TYR'.or.res(I).eq.'PHE'.or.res(I).eq.
1303      1  'TRP').and.aname(I).eq.' CD1')then
1304         iat1=I
1305         do 61 J=-5,5
1306         if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CG ')iat2=I+J
1307         if(ires(I+J).eq.ires(I).and.aname(I+J)(3:4).eq.'E1')iat3=I+J
1308 61      continue
1309         if(iat2.gt.0.and.iat3.gt.0) goto 62
1310        else if((res(I).eq.'TYR'.or.res(I).eq.'PHE'.or.res(I).eq.
1311      1  'HIS').and.aname(I).eq.' CD2') then
1312         iat1=I
1313         do 63 J=-5,5
1314         if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CG ')iat2=I+J
1315         if(ires(I+J).eq.ires(I).and.aname(I+J)(3:4).eq.'E2')iat3=I+J
1316 63      continue
1317         if(iat2.gt.0.and.iat3.gt.0) goto 62
1318        else if((res(I).eq.'TYR'.or.res(I).eq.'PHE').and.aname(I).
1319      1  eq.' CE1') then
1320         iat1=I
1321         do 64 J=-5,5
1322         if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CD1')iat2=I+J
1323         if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CZ ')iat3=I+J
1324 64      continue
1325         if(iat2.gt.0.and.iat3.gt.0) goto 62
1326        else if((res(I).eq.'TYR'.or.res(I).eq.'PHE').and.aname(I).
1327      1  eq.' CE2') then
1328         iat1=I
1329         do 65 J=-5,5
1330         if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CD2')iat2=I+J
1331         if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CZ ')iat3=I+J
1332 65      continue
1333         if(iat2.gt.0.and.iat3.gt.0) goto 62
1334        else if((res(I).eq.'PHE').and.aname(I).eq.' CZ ') then
1335         iat1=I
1336         do 66 J=-5,5
1337         if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CE1')iat2=I+J
1338         if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CE2')iat3=I+J
1339 66      continue
1340         if(iat2.gt.0.and.iat3.gt.0) goto 62
1341        else if(res(I).eq.'HIS'.and.aname(I).eq.' CE1') then
1342         iat1=I
1343         do 67 J=-5,5
1344         if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' ND1')iat2=I+J
1345         if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' NE2')iat3=I+J
1346 67      continue
1347         if(iat2.gt.0.and.iat3.gt.0) goto 62
1348        else if(res(I).eq.'TRP'.and.aname(I).eq.' CE3') then
1349         iat1=I
1350         do 68 J=-8,8
1351         if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CD2')iat2=I+J
1352         if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CZ3')iat3=I+J
1353 68      continue
1354         if(iat2.gt.0.and.iat3.gt.0) goto 62
1355        else if(res(I).eq.'TRP'.and.aname(I).eq.' CZ3') then
1356         iat1=I
1357         do 69 J=-8,8
1358         if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CE3')iat2=I+J
1359         if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CH2')iat3=I+J
1360 69      continue
1361         if(iat2.gt.0.and.iat3.gt.0) goto 62
1362        else if(res(I).eq.'TRP'.and.aname(I).eq.' CH2') then
1363         iat1=I
1364         do 70 J=-8,8
1365         if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CZ3')iat2=I+J
1366         if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CZ2')iat3=I+J
1367 70      continue
1368         if(iat2.gt.0.and.iat3.gt.0) goto 62
1369        else if(res(I).eq.'TRP'.and.aname(I).eq.' CZ2') then
1370         iat1=I
1371         do 71 J=-8,8
1372         if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CH2')iat2=I+J
1373         if(ires(I+J).eq.ires(I).and.aname(I+J).eq.' CE2')iat3=I+J
1374 71      continue
1375         if(iat2.gt.0.and.iat3.gt.0) then
1376 62      iproton=iproton+1
1377         call addsp2(x,y,z,xh,yh,zh,iproton,iat1,iat2,iat3,2)
1378         anameh(iproton)=' H  '
1379         anameh(iproton)(3:4)=aname(iat1)(3:4)
1380         resh(iproton)=res(iat1)
1381         iresh(iproton)=ires(iat1)
1382         endif
1383        else
1384 C presumably a heavy atom that doesn't have protons, or one
1385 C with exchangeable protons that I can't be bothered to do
1386        endif
1387 10     continue
1388        return
1389        end
1390 C===================================================================
1391         subroutine addsp2(x,y,z,xh,yh,zh,Ih,Iat1,Iat2,Iat3,Ibond)
1392         DIMENSION x(5000),y(5000),z(5000)
1393         DIMENSION xh(3000),yh(3000),zh(3000)
1394         DIMENSION r12(3),u12(3),r13(3),u13(3),r14(3),u14(3)
1395 C If Ibond=1, it's a NH bond, if Ibond=2 it's a CH (aromatic) bond
1396         if(Ibond.eq.1)blength=0.98
1397         if(Ibond.eq.2)blength=1.08
1398         r12(1)=x(Iat1)-x(Iat2)
1399         r12(2)=y(Iat1)-y(Iat2)
1400         r12(3)=z(Iat1)-z(Iat2)
1401         call unitv(r12,u12,d)
1402         r13(1)=x(Iat1)-x(Iat3)
1403         r13(2)=y(Iat1)-y(Iat3)
1404         r13(3)=z(Iat1)-z(Iat3)
1405         call unitv(r13,u13,d)
1406         r14(1)=u12(1)+u13(1)
1407         r14(2)=u12(2)+u13(2)
1408         r14(3)=u12(3)+u13(3)
1409         call unitv(r14,u14,d)
1410         xh(Ih)=x(Iat1)+blength*u14(1)
1411         yh(Ih)=y(Iat1)+blength*u14(2)
1412         zh(Ih)=z(Iat1)+blength*u14(3)
1413         return
1414         end
1415 C------------------------------------------------------------------
1416         subroutine addonesp3(x,y,z,xh,yh,zh,Ih,Iat1,Iat2,Iat3,Iat4)
1417         DIMENSION x(5000),y(5000),z(5000)
1418         DIMENSION xh(3000),yh(3000),zh(3000)
1419         DIMENSION r12(3),u12(3),r13(3),u13(3),r14(3),u14(3)
1420         DIMENSION r15(3),u15(3)
1421         blength=1.08
1422         r12(1)=x(Iat1)-x(Iat2)
1423         r12(2)=y(Iat1)-y(Iat2)
1424         r12(3)=z(Iat1)-z(Iat2)
1425         call unitv(r12,u12,d)
1426         r13(1)=x(Iat1)-x(Iat3)
1427         r13(2)=y(Iat1)-y(Iat3)
1428         r13(3)=z(Iat1)-z(Iat3)
1429         call unitv(r13,u13,d)
1430         r14(1)=x(Iat1)-x(Iat4)
1431         r14(2)=y(Iat1)-y(Iat4)
1432         r14(3)=z(Iat1)-z(Iat4)
1433         call unitv(r14,u14,d)
1434         r15(1)=u12(1)+u13(1)+u14(1)
1435         r15(2)=u12(2)+u13(2)+u14(2)
1436         r15(3)=u12(3)+u13(3)+u14(3)
1437         call unitv(r15,u15,d)
1438         xh(Ih)=x(Iat1)+blength*u15(1)
1439         yh(Ih)=y(Iat1)+blength*u15(2)
1440         zh(Ih)=z(Iat1)+blength*u15(3)
1441         return
1442         end
1443 C------------------------------------------------------------------
1444         subroutine addtwosp3(x,y,z,xh,yh,zh,Ih,Iat1,Iat2,Iat3)
1445 C This is based on the GEN routine in VNMR (J. Hoch)
1446         DIMENSION x(5000),y(5000),z(5000)
1447         DIMENSION xh(3000),yh(3000),zh(3000)
1448         DIMENSION r12(3),r13(3),side(3),us(3)
1449         DIMENSION add(3),uadd(3),sub(3),usub(3)
1450         blength=1.08
1451         alpha=54.75*3.14159/180.   !H-C-H=109.5
1452         r12(1)=x(Iat1)-x(Iat2)
1453         r12(2)=y(Iat1)-y(Iat2)
1454         r12(3)=z(Iat1)-z(Iat2)
1455         r13(1)=x(Iat1)-x(Iat3)
1456         r13(2)=y(Iat1)-y(Iat3)
1457         r13(3)=z(Iat1)-z(Iat3)
1458         sub(1)=r12(1)-r13(1)
1459         sub(2)=r12(2)-r13(2)
1460         sub(3)=r12(3)-r13(3)
1461         call unitv(sub,usub,d)
1462         add(1)=r12(1)+r13(1)
1463         add(2)=r12(2)+r13(2)
1464         add(3)=r12(3)+r13(3)
1465         call unitv(add,uadd,d)
1466         call vprod(usub,uadd,side,sintemp,tmp)
1467         call unitv(side,us,d)
1468         xh(Ih-1)=x(Iat1)+blength*(uadd(1)*cos(alpha)-us(1)*sin(alpha))
1469         yh(Ih-1)=y(Iat1)+blength*(uadd(2)*cos(alpha)-us(2)*sin(alpha))
1470         zh(Ih-1)=z(Iat1)+blength*(uadd(3)*cos(alpha)-us(3)*sin(alpha))
1471         xh(Ih)=x(Iat1)+blength*(uadd(1)*cos(alpha)+us(1)*sin(alpha))
1472         yh(Ih)=y(Iat1)+blength*(uadd(2)*cos(alpha)+us(2)*sin(alpha))
1473         zh(Ih)=z(Iat1)+blength*(uadd(3)*cos(alpha)+us(3)*sin(alpha))
1474         return
1475         end
1476 C------------------------------------------------------------------
1477         subroutine addthreesp3(x,y,z,xh,yh,zh,Ih,Iat1,Iat2,Iat3)
1478 C This is based on the GEN routine in VNMR (J. Hoch)
1479         DIMENSION x(5000),y(5000),z(5000)
1480         DIMENSION xh(3000),yh(3000),zh(3000),twist(3),utw(3)
1481         DIMENSION r12(3),u12(3),r23(3),perp(3),side(3),uside(3)
1482         REAL perp
1483         blength=1.08
1484         beta=70.5*3.14159/180.  !180-109.5
1485         cosbeta=cos(beta)
1486         sinbeta=sin(beta)
1487         cosgam=0.866     !cos(30)
1488 c The methyl protons should be staggered to Iat3
1489         r12(1)=x(Iat1)-x(Iat2)
1490         r12(2)=y(Iat1)-y(Iat2)
1491         r12(3)=z(Iat1)-z(Iat2)
1492         call unitv(r12,u12,d)
1493         r23(1)=x(Iat2)-x(Iat3)
1494         r23(2)=y(Iat2)-y(Iat3)
1495         r23(3)=z(Iat2)-z(Iat3)
1496         vert=(r12(1)*r23(1)+r12(2)*r23(2)+r12(3)*r23(3))/d
1497         perp(1)=u12(1)*vert
1498         perp(2)=u12(2)*vert
1499         perp(3)=u12(3)*vert
1500         side(1)=r23(1)-perp(1)
1501         side(2)=r23(2)-perp(2)
1502         side(3)=r23(3)-perp(3)
1503         call unitv(side,uside,d)
1504         call vprod(u12,uside,twist,sintemp,tmp)
1505         call unitv(twist,utw,d)
1506         u12(1)=u12(1)*cosbeta
1507         u12(2)=u12(2)*cosbeta
1508         u12(3)=u12(3)*cosbeta
1509         xh(Ih-2)=x(Iat1)+blength*(u12(1)+uside(1)*sinbeta)
1510         yh(Ih-2)=y(Iat1)+blength*(u12(2)+uside(2)*sinbeta)
1511         zh(Ih-2)=z(Iat1)+blength*(u12(3)+uside(3)*sinbeta)
1512         xh(Ih-1)=x(Iat1)+blength*(u12(1)-sinbeta*(
1513      &   uside(1)/2.0-utw(1)*cosgam))
1514         yh(Ih-1)=y(Iat1)+blength*(u12(2)-sinbeta*(
1515      &   uside(2)/2.0-utw(2)*cosgam))
1516         zh(Ih-1)=z(Iat1)+blength*(u12(3)-sinbeta*(
1517      &   uside(3)/2.0-utw(3)*cosgam))
1518 C NB The /2.0 is actually *sin(30)
1519         xh(Ih)=x(Iat1)+blength*(u12(1)-sinbeta*(
1520      &   uside(1)/2.0+utw(1)*cosgam))
1521         yh(Ih)=y(Iat1)+blength*(u12(2)-sinbeta*(
1522      &   uside(2)/2.0+utw(2)*cosgam))
1523         zh(Ih)=z(Iat1)+blength*(u12(3)-sinbeta*(
1524      &   uside(3)/2.0+utw(3)*cosgam))
1525         return
1526         end