src/contrib/total.f

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