c Read in coordinates and calculate solvent accessible surface area
c 
c To calculate area, this program will distribute a given number of points evenly on the surface of each atom.
c Points will be distributed on the path of a spherical spiral, from north to south pole (-Z to +Z).
c Write point coordinates as atoms in .pdb file format.
c Using Golden Spiral algorithm described by Christopher Carlson.
c http://blog.wolfram.com/2011/07/28/how-i-made-wine-glasses-from-sunflowers/
c Coding in python described here:
c https://www.mcadcentral.com/threads/25421-Even-Distribution-of-Points-on-a-Sphere
c Accracy described in Fig 24 of the following reference
c https://scholar.rose-hulman.edu/rhumj/vol18/iss2/5/
c With this Golden algorithm,aAccuracy increases asymptotically with N. 
c When N=800, accuracy is better than 84%.
c This golden spiral is preferable to the Jonathan Kogan algorithm where accuracy varies wildly with N.
c Program may give spurious results if you have a resid zero. Please renumber.

c Enter pH. Calculate solvation energy.
c Atomic Solvation Parameters (ASP)
c Eisenberg D, Wesson M, Yamashita M (1989) 
c Interpretation of protein folding and binding with atomic solvation parameters. 
c Chemica Scripta 29A: 217–221.
c Terminal nitrogens and oxygens should be given residue names NTE and CTE to take into account
c the possibility of charge. See below. 
c
c The B-factors in the output pdb file will contain the stabilization energies in kcal/mol. 
c Have fun! Mikey
c
c Note: Bonded residues within a chain must be conseculatively numbered. 
c For example, spliced tau in Picks disease must be renumbered so residues are consecutive. 
c No gaps in residue range, unless the residues are truly not bonded to each other.
c  Michael R. Sawaya June 6, 2018
c
c Version 02  kluge KLUGE kluge KLUGE kluge KLUGE kluge KLUGE kluge
c     Calculate solvation free energy, using assumption all chargable atoms
c     are charged if SASAfold.gt.thresh, and use a kluge to reverse sign 
c     of the atomic solvation parameter (A.S.P.) if the culprit of the atom/ion 
c     is a complementary charge. Apply reversed A.S.P. only to the area buried 
c     in the atom-culprit contact. To all other culprits of the atom, use the
c     uncharged A.S.P. regardless of how much SASA is exposed on the atom, as
c     long as the oppositely charged atoms are within 5A distance of each other.
c     Version 01 did not have the SASA override based on 5 A proximity.
c     Therefore version 02 gives more favorable stabilization energy.
c
c Version 03 Add a routine to set A.S.P = 0 for buried polar atoms in 
c     Asn, Gln, Asp, Glu, Ser, His, Lys, Arg, Thr, Trp, Tyr if they are within H-bond distance
c     of another heteroatom.
c
c
c Version 03.1 Can handle charged N and C-termini if resn is explicitly labeled NTE or CTE in pdb file as follows.
c ATOM      1  N   NTE A   1      18.009   4.003   4.619  1.00  6.11           N
c ATOM     88  O   CTE A   6       1.748   1.037   3.962  1.00  7.90           O
c ATOM     90  OXT CTE A   6       0.234   1.134   5.458  1.00  9.06           O
c and Cterm oxygen must be labeled OXT as shown above (not OT)
c Mikey Oct 23, 2018
c
c Version 03.2 Prints areaculprita to log instead of fxnarea. Minor point. 
c Helpful because areaculprita is currently used to calculate kluge energy not fxnarea.
c
c Version 03.3 is like 03.2 but re-written kluge code to be shorter.
c
c Version 03.4 imposes distance cutoff of 3.5 A in order for ion pair to be assigned asp=+37/+38. 
c All other culprits, including ion pairs with distance greater than 3.5A are assigned asp=0.
c In previous version all ion pairs between an atom and its culprit were assigned asp=+37/+38
c This version gives lower stabilization energies. It is contrary to Guenther NFGAFS expts.
c
c Version 4.0 Uses standard A.S.P.s except for the following cases:
c 1) A.S.P. = 0 for sidechain N and O of Asn or Gln that have two hydrogen bonds (otherwise use std negative A.S.P)
c 2) A.S.P. = standard negative values for sidechain N and O of non-Asn/Gln residues. No exceptions.
c 3) A.S.P. = -9 for charged N or O involved in ion pair.
c 4) A.S.P. = -37/-38 for charged O or N not involved in ion pair (distance between charged atoms < 3.5 A)
c 5) A.S.P. = 0 for backbone atoms in hydrogen bond.
c Algorithm suggested by DSE on 2018-11-13 
c Also assigns charge state explicitly to individual atoms of binuclear ions like Asp, Glu, and trinuclear Arg.
c
c Version 5.0 imposes additional conditions to merit the energetic bonuses of Version 4.0. Ammendments are as follows:
c 1) A.S.P. = 0 for sidechain N and O of Asn or Gln that have two hydrogen bonds ***AND*** 
c    area buried greater than 40 A^2 (otherwise use std negative A.S.P)
c    That is, a polar ladder is insufficient to merit the bonus of A.S.P=0; the Asn/Gln N and O must also be buried.
c 3) A.S.P. = -9 for charged N or O involved in ion pair ***AND*** the atoms have area buried greater than 50 A^2.
c    That is, an ion pair is insufficient to merit the bonus of A.S.P=-9; the ionized atoms must also be buried.
c    50 A^2 is very high area buried requirement to ensure the ion pair is not on surface or near solvent channel.
c Algorithm suggested by DSE on 2018-12-03 
c I changed distance cutoff for ion pairs from 3.5 to 3.7 to allow for errors in coordinates. e.g. 6gx5 Asp348-Lys298 pair.
c
c Version 5.1 has more reliable indicator of when Asp,Glu,Lys,Arg,Asn,Gln are buried. Instead of areaburied >50 or 40 A^2,
c use threshhold of sasafolded < 5.0 A^2 for ions, and 5 A^2 for Gln/Asn.
c
c Version 5.2 increases the threshhold SASAfolded for ion pairs from 5 A^2 to 40 A^2
c in order to give energetic bonus to R298-E315 pair on TDP-43 (PDB ID 6n3c) which was rather solvent exposed on one side. 
c Also, change distance cutoff for ion pairs from 3.7 to 4.2 A to qualify as ion pair.
c
c Version 5.3 counts residues and gives energy per residue. 2018-12-31
c
c Version 5.4 sounds an alarm if atoms are missing from any side chain. 2019-01-02
c
c Version 6.0 includes an entropy penalty. Outputs pdb with sfe+entropy in B-factor column. 
c Uses Pickett and Sternberg 1993 sidechain entropy terms
c
c Version 6.1 Replace entropy terms from Pickett and Sternberg 1993 with 
c mean values from 6 models. Doig and Sternberg 1995.
c
c Version 6.2 Scale TdS from Doig and Sternberg 1995 by a factor of 0.704 to make sfe+TdS for Fus = 0
c
c Version 6.3 Scale TdS from Doig and Sternberg 1995 by percent area buried.
c percent area buried = (sasaisol(i)-sasafold(i))/sasaisol(i) for each atom i.
c
c Version 6.4 Scale TdS from Doig and Sternberg 1995 by 1.0 if SASAfolded.gt.5 A^2
c Scale by 0.850 for atoms if SASAfolded.ge.5 A^2
c
c Version 6.5 Use TdS values from Koehl and Delarue 1994 because DSE thinks these are more reasonable than Doig and Sternberg 1995.
c DSE objects to large penalty to SER residues incurred by Sternberg and abated by Delarue.
c
c Version 6.7 Scale TdS values from Koehl and Delarue 1994 by the percentage of side chain surface area buried in assembly.
c Previously, entropic contributions for each amino acid type were given the same value regardless of whether the side chain 
c is buried or exposed. For example, buried lysine side chains incur the same penalty as exposed lysine side chains. 
c This uniform treatment of amino acid types seems to be an oversimplification. We expect that a wider range of conformations 
c are available to  surface exposed side chains compared to buried side chains. Therefore, solvent exposed side chains should 
c incur less entropic penalty than buried side chains.
c
c Version 6.8 Scale TdS values from Koehl and Delarue 1994 by 1.2 * (percentage of side chain surface area buried in assembly)^2.
c Further lowers penalty on outward facing residues and raises penalty on inward facing residues.
c
c Version 7.0 Give more enthalpic stabilization to ion pairs in the following two ways:
c 1) Raise threshold SASA folded from 40 A^2 (kluge5.2) to 50 A^2 (to benefit ion pair on tau PHF 5o3l K343-D345).
c 2) Boost energetic bonus for ion pairs from A.S.P.=-9 to A.S.P.=0 depending on distance between acceptor and donor.
c    A.S.P.=-9* ((dist-2.8A)/2.8A)^2.  
c    Since ion pairs are considered only when dist < 4.2A (kluge5.2) the range is -4.5 < A.S.P .<0
c
c Version 7.1 Raise threshold for ion pair distance from 4.2A to 4.5A. It made 5o3l and 5o3t closer in energy as expected.
c
c Version 7.2a go back to simpler TDs scaling from version 6.7. Koehl&Delarue1994 * %sidechainsurfaceareaburied. 
c It made 5o3l and 5o3t even closer in energy, as expected. It punished 6qjh heparin induced tau snake polymorph a little, as desired.
c
c Version 7.2c Modify to include uncharged amino NH2 on amidated IAPP. Atom name should be NXT, and residue name should be ATE (amide terminus).
c ATOM   2968  NXT ATE A  37      87.371  81.623  93.423  1.00 41.81           N
c
c Version 7.2d Modify threshhold distance of Hbond from <3.2 to <3.3 A. Line 1231
c
       character pdbinfile1*100
c      Output PDB file containing point coordinates distributed evenly on surface of a sphere
       character dotfile*100 
c      lineinput is a temporary storage of pdb input lines to discriminate atoms from header
       character lineinput*80
c      Chain id of contacting molecules
       character ch(9000)*1
c      Resid of contacting molecules
       integer resid(9000)
c      Number of atoms read from PDB file.
       integer n
c      Counting how many chains there are in input PDB.
       character prevchain*1,chainlist(62)*1
       integer nchains
c      PDB items to read and parse.
       real x(9000),y(9000),z(9000)
       character resn(9000)*3,atnam(9000)*5
c      Number of points to be equally distributed on the surface of a sphere
       integer Npoints
c      Various parameters for spiral calculation
       real r, off, yy, inc
c      Angle in radians for spherical spiral calculation.
       real phi
c      Coordinates of points distributed on surface of a unit sphere.
       real xC(1000),yC(1000),zC(1000)
c      Coordinates of points on surface of each atom
       real xdot(9000,1000),ydot(9000,1000),zdot(9000,1000)
c      Record of each resid that shields a dot from solvent access.
       integer residculpritd(9000,1000,20)
c      Record of each resname that shields solvent from a given dot.
       character*3 resnculpritd(9000,1000,20) 
       character*1 chculpritd(9000,1000,20)
       character*3 atnamculpritd(9000,1000,20)
c      Number of atoms that encompass a dot.
       integer nculpritd(9000,1000)
c      Solvent accessible area of an atom in isolated tripeptide
       real sasaisol(9000)
c      Solvent accessible area of an atom in context of fold/complex.
       real sasafold(9000)
c      Record of each resid that shields an atom from solvent access.
       integer residculprita(9000,85)
       integer pct25, pct50, pct75
c      Record of each resname that shields solvent from a given dot.
       character*3 resnculprita(9000,85)
       character*1 chculprita(9000,85)
       character*3 atnamculprita(9000,85)
c      Number of atoms that tresspass on an atom.
       integer nculprita(9000)
c      Area of an atom that is shielded by each of up to 20 culprits.
       real areaculprita(9000,85)
c      Solvent accessible area difference between isolated and folded state.
       real aburied(9000)
c      Solvation Free Energy total
       real sumnrgy
c      pH for solvation Free energy calculation. 
       real pHvalu
c      Solvation Free Energy for an atom
       real sfe(9000)
c      Sum the area of the victim buried by all its culprits.
c      Will be a bigger number than actual area buried because 
c      one dot of victim surface area can be enclosed by multiple atoms.
       real sumareaculprita(9000)
c      Is ionizable atom neutralized by an opposite charge less than 5.0A away?
c      charge(i) = 1 if no opposite charge, charge(i) = 0 if opposite charge is close by.
       real charge(9000)
       real asp(9000)
       character*2  atomtyp, culptyp
       integer nhbond(9000)
c      Number of residues in chain A->Z=26 a->z=26 0->9=10 23+23+10=62
       integer nresid(62)
c      Residue name of previous amino acid
       character prevresn*3
c      Number of each of 20 amino acids in chain j
       integer nuA(62),nuC(62),nuD(62),nuE(62),nuF(62),nuG(62)
       integer nuH(62),nuI(62),nuK(62),nuL(62),nuM(62),nuN(62)
       integer nuP(62),nuQ(62),nuR(62),nuS(62),nuT(62),nuV(62)
       integer nuW(62),nuY(62)
c      TdS entropy penalty summed over each chain j
       real entch(62)
c      TdS entropy penalty per atom
       real entpatom(9000)
c      Nth residue encountered (rescnt is a concatenation of residue count)
       integer rescnt
c      resid of nth residue
       integer nthresresid(9000)
c      chain id of nth residue
       character nthresch(9000)*1
c      sasafold and sasaisol for nth residue
       real nthressasafold(9000),nthressasaisol(9000)
c      Perecnt area buried on Nth side chain for Nth residue
       real nthrespctabu(9000)
c      Perecnt area buried on Nth side chain indexed on Ith atom
       real nthatompctabu(9000)
c      deltaG for Nth residue
       real nthresdG(9000)
c      deltaG of Nth residue indexed on Ith atom
       real nthatomdG(9000)
c
        write(*,*) 'Please type name of .pdb file'
        read(*,'(a100)')pdbinfile1
        open (1,file=pdbinfile1, status='old', form='formatted')
        print*, pdbinfile1
        print*, 'PLEASE give pH value for calculation'
        read (*,"(F4.1)")pHvalu
 	open (3,file='dotfile.pdb',status='new',form='formatted')
 	open (4,file='sasaisolated.pdb',status='new',form='formatted')
 	open (5,file='sasafolded.pdb',status='new',form='formatted')
 	open (7,file='aburied.pdb',status='new',form='formatted')
 	open (8,file='sfesasathres20.pdb',status='new',form='formatted')
 	open (9,file='sfekluge7.2.pdb',status='new',form='formatted')
 	open (10,file='hbonds.pml',status='new',form='formatted')
 	open (11,file='aspkluge7.2.pdb',status='new',form='formatted')
 	open (12,file='pctsideburied.pdb',status='new',form='formatted')
 	open (13,file='sfekluge7.2.sumperresid.pdb',status='new',
     1  form='formatted')
 	open (14,file='hbunsatisfied.pdb',status='new',form='formatted')
c       If value of sasafold exceeds thresh, then ionizable atom is charged
        thresh=20
c       Now I am reading in the coordinates from pdbinfile1.
        prevchain='.'
        nchains=0
        nresid(1)=0
        n1statom=0
        natomprevresid=0
5       n=1
2       format(a80)
4       read (1,2,end=6)lineinput
        if((lineinput(1:4).ne.'ATOM').and.(n1statom.eq.0))then
         write(3,2)lineinput
         goto 4
        endif
c       skip over END and TER records.
        if((lineinput(1:4).ne.'ATOM').and.(n1statom.ne.0))goto 4
        n1statom=1
        if(n.gt.9000)then
         print*, 'Redimension for more than 9000 atoms.'
         go to 999
        endif
        atnam(n)=lineinput(13:17)  
        resn(n)=lineinput(18:20)
        ch(n)=lineinput(22:22)
        read(lineinput(23:26),'(i4)')resid(n)
        read(lineinput(31:38),*)x(n)
        read(lineinput(39:46),*)y(n)
        read(lineinput(47:54),*)z(n)
        if(n.eq.1)goto 312
        if((ch(n).eq.prevchain).and.(resid(n).eq.resid(n-1)))then
         natomprevresid=natomprevresid+1
         if(atnam(n).eq.' CA  ')prevresn=resn(n)
         if(atnam(n).eq.' OXT ')natomprevresid=natomprevresid-1
         if(atnam(n).eq.' NXT ')natomprevresid=natomprevresid-1
        endif
        if((ch(n).eq.prevchain).and.(resid(n).ne.resid(n-1)))then
         nresid(nchains)=nresid(nchains)+1
c        print*,ch(n),resid(n-1),prevresn,natomprevresid
         if((prevresn.eq.'GLY').and.(natomprevresid.ne.4))goto 998
         if((prevresn.eq.'ALA').and.(natomprevresid.ne.5))goto 998
         if((prevresn.eq.'SER').and.(natomprevresid.ne.6))goto 998
         if((prevresn.eq.'CYS').and.(natomprevresid.ne.6))goto 998
         if((prevresn.eq.'VAL').and.(natomprevresid.ne.7))goto 998
         if((prevresn.eq.'THR').and.(natomprevresid.ne.7))goto 998
         if((prevresn.eq.'PRO').and.(natomprevresid.ne.7))goto 998
         if((prevresn.eq.'ILE').and.(natomprevresid.ne.8))goto 998
         if((prevresn.eq.'LEU').and.(natomprevresid.ne.8))goto 998
         if((prevresn.eq.'ASP').and.(natomprevresid.ne.8))goto 998
         if((prevresn.eq.'ASN').and.(natomprevresid.ne.8))goto 998
         if((prevresn.eq.'MET').and.(natomprevresid.ne.8))goto 998
         if((prevresn.eq.'GLU').and.(natomprevresid.ne.9))goto 998
         if((prevresn.eq.'GLN').and.(natomprevresid.ne.9))goto 998
         if((prevresn.eq.'LYS').and.(natomprevresid.ne.9))goto 998
         if((prevresn.eq.'HIS').and.(natomprevresid.ne.10))goto 998
         if((prevresn.eq.'PHE').and.(natomprevresid.ne.11))goto 998
         if((prevresn.eq.'ARG').and.(natomprevresid.ne.11))goto 998
         if((prevresn.eq.'TYR').and.(natomprevresid.ne.12))goto 998
         if((prevresn.eq.'TRP').and.(natomprevresid.ne.14))goto 998
         natomprevresid=1
         if(atnam(n).eq.' CA  ')prevresn=resn(n)
        endif
312     if(ch(n).ne.prevchain)then
         nchains=nchains+1
         if(nchains.gt.62)then
          print*, 'Redimension for more than 62 chains.'
          goto 999
         endif
         nresid(nchains)=1
         chainlist(nchains)=ch(n)
         prevchain=ch(n)
c        if(n.ne.1)print*,ch(n-1),resid(n-1),prevresn,natomprevresid
         if((prevresn.eq.'GLY').and.(natomprevresid.ne.4))goto 998
         if((prevresn.eq.'ALA').and.(natomprevresid.ne.5))goto 998
         if((prevresn.eq.'SER').and.(natomprevresid.ne.6))goto 998
         if((prevresn.eq.'CYS').and.(natomprevresid.ne.6))goto 998
         if((prevresn.eq.'VAL').and.(natomprevresid.ne.7))goto 998
         if((prevresn.eq.'THR').and.(natomprevresid.ne.7))goto 998
         if((prevresn.eq.'PRO').and.(natomprevresid.ne.7))goto 998
         if((prevresn.eq.'ILE').and.(natomprevresid.ne.8))goto 998
         if((prevresn.eq.'LEU').and.(natomprevresid.ne.8))goto 998
         if((prevresn.eq.'ASP').and.(natomprevresid.ne.8))goto 998
         if((prevresn.eq.'ASN').and.(natomprevresid.ne.8))goto 998
         if((prevresn.eq.'MET').and.(natomprevresid.ne.8))goto 998
         if((prevresn.eq.'GLU').and.(natomprevresid.ne.9))goto 998
         if((prevresn.eq.'GLN').and.(natomprevresid.ne.9))goto 998
         if((prevresn.eq.'LYS').and.(natomprevresid.ne.9))goto 998
         if((prevresn.eq.'HIS').and.(natomprevresid.ne.10))goto 998
         if((prevresn.eq.'PHE').and.(natomprevresid.ne.11))goto 998
         if((prevresn.eq.'ARG').and.(natomprevresid.ne.11))goto 998
         if((prevresn.eq.'TYR').and.(natomprevresid.ne.12))goto 998
         if((prevresn.eq.'TRP').and.(natomprevresid.ne.14))goto 998
         natomprevresid=1
         if(atnam(n).eq.' N   ')prevresn=resn(n)
        endif
        n=n+1
        go to 4
6       n=n-1
        print*, n, 'atoms read from input PDB file.'
        print*, nchains, 'chains found.'
        write(6,16)chainlist(1),chainlist(nchains) 
16      format('First chain=',a1,'. Last chain=',a1,'.')
        do 313, icha=1,nchains
         write(6,314)nresid(icha),chainlist(icha)
314      format(i5,' residues in chain ',a1)
313     continue
c
c      Calculate coordinates of points on a unit sphere.
c      Use 1000 points for all elements. Plan to adjust radii in paragraph below.
       print*,'Start calculating 1000 points on surface of unit sphere.'
       Npoints=1000
       inc=3.14159 *(3.0-sqrt(5.0))
       off=2.0/Npoints
       do 101, k=1,Npoints
        yy=((float(k-1))*off)-1.0+(off/2.0)
        r=sqrt(1-(yy*yy))
        phi=(k-1)*inc
        xC(k)=r*(cos(phi))
        yC(k)=r*sin(phi)
        zC(k)=yy
c       print*, k-1, yy, r, phi, xC(k),yC(k),zC(k)
101    continue
       print*, 'Done calculating 1000 points on surface of unit sphere.'
       

c      Make dot surface
       print*, 'Start assigning dot surface to each atom.'
       proberadius=1.4
c      van der Waal radii for carbon, nitrogen, oxygen, and sulfur. Values match areaimol.
       radC=1.75
       radN=1.60
       radO=1.55
       radS=1.80
       do 40, i=1,n
        if(atnam(i)(2:2).eq.'C')then
         radatom=radC+proberadius
         goto 37
        endif
        if(atnam(i)(2:2).eq.'N')then
         radatom=radN+proberadius
         goto 37
        endif
        if(atnam(i)(2:2).eq.'O')then
         radatom=radO+proberadius
         goto 37
        endif
        if(atnam(i)(2:2).eq.'S')then
         radatom=radS+proberadius
         goto 37
        endif
        print*, 'No definition for element',atnam(i)
        goto 999
37      do 41, k=1,Npoints
        xdot(i,k)=x(i)+(xC(k)*radatom)
        ydot(i,k)=y(i)+(yC(k)*radatom)
        zdot(i,k)=z(i)+(zC(k)*radatom)
41      continue        
40     continue
       print*, 'Done assigning dot surface to each atom.'

c      Initialize dot records.
       print*, 'Start initializing culprit array for dots.'
       pct25=0.25*n*Npoints*20
       pct50=0.50*n*Npoints*20
       pct75=0.75*n*Npoints*20
       current=0
       do 57, i=1,n
        do 58, k=1,Npoints
          nculpritd(i,k)=0
         do 59, l=1,20
          current=current+1
          if(current.eq.pct25)print*,'25% done initializing'
          if(current.eq.pct50)print*,'50% done initializing'
          if(current.eq.pct75)print*,'75% done initializing'
          residculpritd(i,k,l)=0
          resnculpritd(i,k,l)='non'
          chculpritd(i,k,l)='0'
          atnamculpritd(i,k,l)='non'
59       continue
58      continue
57     continue
       print*, 'Finished initializing culprit array for dots.'
c
c      We want to know which dots on an atom's surface are enveloped by other atoms.
c      Flag dots that interpenetrate atom surfaces, if resid j = resid i, i+1 or i-1  
c      Record which atom is shielding each dot.
       print*, 'Start flagging interpenetrating surface dots.'
       do 45, i=1,n
        do 44, j=1,n
         if(i.eq.j)goto 44
c        Calculate distance between atom i and atom j
         dx=x(i)-x(j)
         dy=y(i)-y(j)
         dz=z(i)-z(j)
         dist=sqrt((dx*dx)+(dy*dy)+(dz*dz))
c        If the distance between atoms is less than the sum of vanderwaal radii, proceed to test for interpenetrating dots.
         if(atnam(i)(2:2).eq.'C')radiusi=radC+proberadius
         if(atnam(i)(2:2).eq.'N')radiusi=radN+proberadius
         if(atnam(i)(2:2).eq.'O')radiusi=radO+proberadius
         if(atnam(i)(2:2).eq.'S')radiusi=radS+proberadius
         if(atnam(j)(2:2).eq.'C')radiusj=radC+proberadius
         if(atnam(j)(2:2).eq.'N')radiusj=radN+proberadius
         if(atnam(j)(2:2).eq.'O')radiusj=radO+proberadius
         if(atnam(j)(2:2).eq.'S')radiusj=radS+proberadius
         sumrad=radiusi+radiusj
         if(dist.lt.sumrad)then
          do 42, k=1,Npoints
c          Testing for interpenetrating dots.
           ddx=xdot(i,k)-x(j)
           ddy=ydot(i,k)-y(j)
           ddz=zdot(i,k)-z(j)
           ddist=sqrt((ddx*ddx)+(ddy*ddy)+(ddz*ddz))
c          If distance from dot(atom i) to nucleus of atom j is shorter than the radius of atom j flag dot.
           if(ddist.lt.radiusj)then
            do 48, l=1,20
             if(residculpritd(i,k,l).eq.0)then
              residculpritd(i,k,l)=resid(j)
              resnculpritd(i,k,l)=resn(j)
              chculpritd(i,k,l)=ch(j)
              atnamculpritd(i,k,l)=atnam(j)(2:4)
              nculpritd(i,k)=nculpritd(i,k)+1
              goto 42
             endif
             if(l.eq.20)then
              print*, 'You got 20 atoms encircling this dot. Crazy!'
               do 83, m=1,20
              print*, ch(i),resid(i),resn(i),atnam(i),chculpritd(i,k,m),
     1    residculpritd(i,k,m),resnculpritd(i,k,m),atnamculpritd(i,k,m)
              print*, 'Redimension variables that end in culpritd'
83             continue
              goto 999
             endif
48          continue
           endif
42        continue
         endif
44      continue
45     continue
       print*, 'Done flagging interpenetrating surface dots.'

c     Tally solvent accessible surface area of "unfolded" chains
      do 47, j=1,nchains
      accessarea=0
      do 50, i=1,n
       if(ch(i).ne.chainlist(j))goto 50
       if(atnam(i)(2:2).eq.'C')radiusi=radC+proberadius
       if(atnam(i)(2:2).eq.'N')radiusi=radN+proberadius
       if(atnam(i)(2:2).eq.'O')radiusi=radO+proberadius
       if(atnam(i)(2:2).eq.'S')radiusi=radS+proberadius
       sasaisol(i)=0.0
       do 51, k=1,Npoints
c       Not all points on an atom's sphere contribute to its surface area.
c       Omit counting area of a point if it is shielded from solvent by culprits in
c       the same residue or adjacent in sequence, as coded below:
c       1) The cuplrit belongs to the same residue as the dot in question.
c          That is, the dot is enclosed in the sphere of another atom (the culprit) in the same residue.
c                              OR
c       2) The culprit is a mainchain N from the next residue.
c                              OR
c       3) The culprit is a mainchain C from the previous residue.
c                              OR
c       4) The dot is from a mainchain atom and the culprit is a mainchain atom from the previous or next residue.
        do 56, l=1,20
c        condition 1
         if((residculpritd(i,k,l).ne.0).and.
     1   (ch(i).eq.chculpritd(i,k,l)).and.
     1   (resid(i).eq.residculpritd(i,k,l)))goto 51
c        condition 2
         if((residculpritd(i,k,l).ne.0).and.
     1   (ch(i).eq.chculpritd(i,k,l)).and.
     1   (resid(i)-residculpritd(i,k,l).eq.1).and.
     1   (atnamculpritd(i,k,l)(1:2).eq.'N '))goto 51
c        condition 3
         if((residculpritd(i,k,l).ne.0).and.
     1   (ch(i).eq.chculpritd(i,k,l)).and.
     1   (resid(i)-residculpritd(i,k,l).eq.-1).and.
     1   (atnamculpritd(i,k,l)(1:2).eq.'C '))goto 51
c        condition 4
         atomtyp='sc'
         if(atnam(i)(2:3).eq.'N ')atomtyp='mc'
         if(atnam(i)(2:3).eq.'C ')atomtyp='mc'
         if(atnam(i)(2:3).eq.'CA')atomtyp='mc'
         if(atnam(i)(2:3).eq.'O ')atomtyp='mc'
         culptyp='sc'
         if(atnamculpritd(i,k,l)(1:2).eq.'N ')culptyp='mc'
         if(atnamculpritd(i,k,l)(1:2).eq.'C ')culptyp='mc'
         if(atnamculpritd(i,k,l)(1:2).eq.'CA')culptyp='mc'
         if(atnamculpritd(i,k,l)(1:2).eq.'O ')culptyp='mc'
c        print*,'atomtyp=',atomtyp,' culptyp=',culptyp,' '
c    1   ,atnam(i)(2:3),' ',
c    1   atnamculpritd(i,k,l)(1:2)
         if((residculpritd(i,k,l).ne.0).and.
     1   (ch(i).eq.chculpritd(i,k,l)).and.
     1   (abs(resid(i)-residculpritd(i,k,l)).eq.1).and.
     1   (atomtyp.eq.'mc').and.(culptyp.eq.'mc'))then
c          print*, 'mainmain',atnam(i)(2:3),atnamculpritd(i,k,l)(1:2)
          goto 51
         endif
56      continue
c       If you made it here, this dot counts as solvent accessible
        accessarea=accessarea+(4.0*3.1415*radiusi*radiusi/Npoints)
        sasaisol(i)=sasaisol(i)+(4.0*3.1415*radiusi*radiusi/Npoints)
51     continue
       write(4,49)i,atnam(i),resn(i),ch(i),resid(i),
     1 x(i),y(i),z(i),sasaisol(i)
49     format('ATOM',i7,1x,a5,a3,1x,a1,i4,4x,3f8.3,'  1.00',f6.2)
50    continue
       write(6,46)chainlist(j),accessarea
46     format('Solvent Accessible Surface Area for chain ',a1,
     1 ' unfolded=' f9.1,' angstroms squared.')
47    continue

c     Intialize more arrays
       print*, 'Start initializing culprit array for atoms.'
       do 67, i=1,n
         nculprita(i)=0
         do 69, m=1,85
          residculprita(i,m)=0
          resnculprita(i,m)='non'
          chculprita(i,m)='0'
          atnamculprita(i,m)='non'
          areaculprita(i,m)=0.0
69       continue
67     continue
       print*, 'Done initializing culprit array for atoms.'

c     Tally area for folded chain.
      do 68, j=1,nchains
      accessarea=0
      do 61, i=1,n
       if(ch(i).ne.chainlist(j))goto 61
c      write(6,73),i,ch(i),resid(i),resn(i),atnam(i)
73     format(i4,1x,a1,1x,i4,1x,a3,1x,a3)
       if(atnam(i)(2:2).eq.'C')radiusi=radC+proberadius
       if(atnam(i)(2:2).eq.'N')radiusi=radN+proberadius
       if(atnam(i)(2:2).eq.'O')radiusi=radO+proberadius
       if(atnam(i)(2:2).eq.'S')radiusi=radS+proberadius
       sasafold(i)=0.0
       do 62, k=1,Npoints
c       Tally area from any dot that escapes the boundary of another atom.
c       This means any atom with nculpritd=0, i.e. no other atom encloses the dot.
        if(nculpritd(i,k).eq.0)then
         accessarea=accessarea+(4.0*3.1415*radiusi*radiusi/Npoints)
         sasafold(i)=sasafold(i)+(4.0*3.1415*radiusi*radiusi/Npoints)
        endif
c       Tally loss of accessible surface area inflicted by each culprit atom.
        if(nculpritd(i,k).ne.0)then
c        print*, nculpritd(i,k)
         do 65, l=1,20
c         Skip over a point if it is shielded from solvent by 
c         an atom less than one residue away in sequence.
          if((residculpritd(i,k,l).ne.0).and.
     1    (ch(i).eq.chculpritd(i,k,l)).and.
     1    ((resid(i).eq.residculpritd(i,k,l)).or.
     1    ((abs(resid(i)-residculpritd(i,k,l)).eq.1).and.
     1    ((atnamculpritd(i,k,l)(1:2).eq.'CA').or.
     1    (atnamculpritd(i,k,l)(1:2).eq.'N ').or.
     1    (atnamculpritd(i,k,l)(1:2).eq.'O ').or.
     1    (atnamculpritd(i,k,l)(1:2).eq.'CB').or.
     1    (atnamculpritd(i,k,l)(1:2).eq.'C ')))))goto 62
c         Skip over this "l" if culprit resid is "0". 
c         i.e. "l"th culprit not found for this atom.
c         Skip to next point k, because no culprit will be found in higher l for this point k.
          if(residculpritd(i,k,l).eq.0)goto 62
          do 66, m=1,85
c          If this atom has no record of this culprit, create record.
           if(residculprita(i,m).eq.0)then
            residculprita(i,m)=residculpritd(i,k,l)
            resnculprita(i,m)=resnculpritd(i,k,l)
            chculprita(i,m)=chculpritd(i,k,l)
            atnamculprita(i,m)=atnamculpritd(i,k,l)
            areaculprita(i,m)=areaculprita(i,m)+
     1      4.0*3.1415*radiusi*radiusi/Npoints
            nculprita(i)=nculprita(i)+1
c           print*,m,chculprita(i,m),residculprita(i,m),
c    1      resnculprita(i,m), atnamculprita(i,m)
            if(m.eq.85)then
             print*,'More than 85 atoms within vdw+radius of this atom.'
             print*,'Increase dimensions of variables ending in ulprita'
             goto 999
            endif
            goto 65
           endif
c          If this atom does have a record of a culprit matching this dot, append area.
           if((residculprita(i,m)).eq.residculpritd(i,k,l).and.
     1     (resnculprita(i,m).eq.resnculpritd(i,k,l)).and.
     1     (chculprita(i,m).eq.chculpritd(i,k,l)).and.
     1     (atnamculprita(i,m).eq.atnamculpritd(i,k,l)))then
            areaculprita(i,m)=areaculprita(i,m)+
     1      4.0*3.1415*radiusi*radiusi/Npoints
            goto 65
           endif
66        continue
65       continue
        endif
62     continue
       write(5,49)i,atnam(i),resn(i),ch(i),resid(i),
     1 x(i),y(i),z(i),sasafold(i)
61    continue
       write(6,74)chainlist(j),accessarea
74     format('Solvent Accessible Surface Area for chain ',a1,
     1 ' in assembly=' f9.1,' angstroms squared.')
68    continue

c     Write out dots in interface between molecules and distant residues in same molecule.
c     Only write 3rd dot for each atom to manage dot file size.
      dotsaver=2
      do 111, i=1,n
       do 112, k=1,Npoints
        if(nculpritd(i,k).eq.0)goto 112
        do 113, l=1,nculpritd(i,k)
c        Skip over a point if it is shielded from solvent by 
c        an atom less than one residue away in sequence.
         if((ch(i).eq.chculpritd(i,k,l)).and.
     1   ((resid(i).eq.residculpritd(i,k,l)).or.
     1   ((abs(resid(i)-residculpritd(i,k,l)).eq.1).and.
     1   ((atnamculpritd(i,k,l)(1:2).eq.'CA').or.
     1   (atnamculpritd(i,k,l)(1:2).eq.'N ').or.
     1   (atnamculpritd(i,k,l)(1:2).eq.'O ').or.
     1   (atnamculpritd(i,k,l)(1:2).eq.'CB').or.
     1   (atnamculpritd(i,k,l)(1:2).eq.'C ')))))goto 112
c        print*,i,k,l,ch(i),resid(i),atnam(i),
c    1 chculpritd(i,k,l),residculpritd(i,k,l),atnamculpritd(i,k,l)(1:2)
113     continue
        dotsaver=dotsaver+1
        if(dotsaver.eq.3)then
         write(3,43)i,atnam(i),resn(i),ch(i),resid(i),
     1   xdot(i,k),ydot(i,k),zdot(i,k)
43       format('ATOM',i7,1x,a5,a3,1x,a1,i4,4x,3f8.3,'  1.00  0.00')
         dotsaver=1
        endif
112    continue
111   continue


c     Sort the culprits by area they enclose on the victim atom.
c     Biggest culprits first.
      do 81, i=1,n
       if(nculprita(i).le.1)goto 81
c      print*,'nculprita(i)=',nculprita(i)
       do 82, mo=2,nculprita(i)
c       print*, 'mo=',mo
        do 85, m=2,nculprita(i)
c       print*, 'm=',m
c       print*, areaculprita(i,1),areaculprita(i,2),areaculprita(i,3)
         if(areaculprita(i,m).gt.areaculprita(i,m-1))then
          oldculpritarea=areaculprita(i,m-1)
          areaculprita(i,m-1)=areaculprita(i,m)
          areaculprita(i,m)=oldculpritarea
         endif
85      continue
82     continue
81    continue 
      
c     Write out delta area and itemize the areas buried by culprits.
      do 70, i=1,n
       aburied(i)=sasaisol(i)-sasafold(i)
       write(7,71)i,atnam(i),resn(i),ch(i),resid(i),
     1 x(i),y(i),z(i),aburied(i),nculprita(i),
     1 chculprita(i,1),residculprita(i,1),
     1 resnculprita(i,1),atnamculprita(i,1),areaculprita(i,1),
     1 chculprita(i,2),residculprita(i,2),
     1 resnculprita(i,2),atnamculprita(i,2),areaculprita(i,2),
     1 chculprita(i,3),residculprita(i,3),
     1 resnculprita(i,3),atnamculprita(i,3),areaculprita(i,3),
     1 chculprita(i,4),residculprita(i,4),
     1 resnculprita(i,4),atnamculprita(i,4),areaculprita(i,4),
     1 chculprita(i,5),residculprita(i,5),
     1 resnculprita(i,5),atnamculprita(i,5),areaculprita(i,5),
     1 chculprita(i,6),residculprita(i,6),
     1 resnculprita(i,6),atnamculprita(i,6),areaculprita(i,6),
     1 chculprita(i,7),residculprita(i,7),
     1 resnculprita(i,7),atnamculprita(i,7),areaculprita(i,7),
     1 chculprita(i,8),residculprita(i,8),
     1 resnculprita(i,8),atnamculprita(i,8),areaculprita(i,8),
     1 chculprita(i,9),residculprita(i,9),
     1 resnculprita(i,9),atnamculprita(i,9),areaculprita(i,9)
71     format('ATOM',i7,1x,a5,a3,1x,a1,i4,4x,3f8.3,'  1.00',f6.2,i9,
     1 2x,a1,i4,1x,a3,1x,a3,f6.2, 
     1 2x,a1,i4,1x,a3,1x,a3,f6.2, 
     1 2x,a1,i4,1x,a3,1x,a3,f6.2, 
     1 2x,a1,i4,1x,a3,1x,a3,f6.2, 
     1 2x,a1,i4,1x,a3,1x,a3,f6.2, 
     1 2x,a1,i4,1x,a3,1x,a3,f6.2, 
     1 2x,a1,i4,1x,a3,1x,a3,f6.2, 
     1 2x,a1,i4,1x,a3,1x,a3,f6.2, 
     1 2x,a1,i4,1x,a3,1x,a3,f6.2) 
70    continue

c     Sum the area of the victim buried by all its culprits.
c     Will be a bigger number than abudif(i).
c     Need it to scale the area buried by charge pair later.
      do 200, i=1,n
       sumareaculprita(i)=0
       do 201, m=1,nculprita(i)
        sumareaculprita(i)=sumareaculprita(i)+areaculprita(i,m)
201    continue
200   continue

c     Calculate solvation free energy, using mikeys original protocol
c     where atoms are not charged if accessible surface area  is less than 20 A^2. 
      do 93, j=1,nchains
       sumnrgy=0.0
       do 91, i=1,n
        if(ch(i).ne.chainlist(j))goto 91
        aburied(i)=sasaisol(i)-sasafold(i)
        if(atnam(i)(2:2).eq.'C')sfe(i)=aburied(i)*18.0
        if(atnam(i)(2:2).eq.'N')sfe(i)=aburied(i)*(-9.0)
        if(atnam(i)(2:2).eq.'O')sfe(i)=aburied(i)*(-9.0)
        if(atnam(i)(2:2).eq.'S')sfe(i)=aburied(i)*(-5.0)
c       print*, 1,i, atnam(i),atnam(i)(2:2),aburied(i),sfe(i)
c       check for charged atoms
c       posibilities for charged atoms are:
c       His ND1 or NE2  pKa 6.08
c       Tyr OH          pKa 10.1
c       Asp OD1 or OD2  pKa  3.9
c       Glu OE1 or OE2  pKa  4.3
c       Lys NZ          pKa 10.5
c       Arg NH1 or NH2  pKa 12.0
c       histidine
        if((pHvalu.le.6.08).and.(sasafold(i).gt.thresh))then
         if((resn(i).eq.'HIS').and.(atnam(i)(2:4).eq.'ND1'))
     1   sfe(i)=(aburied(i)*(-38.)*0.5)+(aburied(i)*(-9.0)*0.5)
         if((resn(i).eq.'HIS').and.(atnam(i)(2:4).eq.'NE2'))
     1   sfe(i)=(aburied(i)*(-38.)*0.5)+(aburied(i)*(-9.0)*0.5)
        endif
c       lysine
        if((pHvalu.le.10.5).and.(sasafold(i).gt.thresh))then
         if(atnam(i)(2:4).eq.'NZ ')sfe(i)=aburied(i)*(-38.)
        endif
c       arginine
        if((pHvalu.le.10.5).and.(sasafold(i).gt.thresh))then
c        if(atnam(i)(2:3).eq.'NH')print*,'I found a charged Arg'
         if(atnam(i)(2:4).eq.'NH1')sfe(i)=aburied(i)*(-38.)*0.5
         if(atnam(i)(2:4).eq.'NH2')sfe(i)=aburied(i)*(-38.)*0.5
        endif
c       aspartate
        if((pHvalu.ge.3.9).and.(sasafold(i).gt.thresh))then
         if((resn(i).eq.'ASP').and.(atnam(i)(2:4).eq.'OD1'))
     1   sfe(i)=(aburied(i)*(-37.)*0.5)+(aburied(i)*(-9.0)*0.5)
         if((resn(i).eq.'ASP').and.(atnam(i)(2:4).eq.'OD2'))
     1  sfe(i)=(aburied(i)*(-37.)*0.5)+(aburied(i)*(-9.0)*0.5)
        endif
c       glutamate
        if((pHvalu.ge.3.9).and.(sasafold(i).gt.thresh))then
c        if(atnam(i)(2:3).eq.'OE')print*,'I found a charged Glu'
         if((resn(i).eq.'GLU').and.(atnam(i)(2:4).eq.'OE1'))
     1   sfe(i)=(aburied(i)*(-37.)*0.5)+(aburied(i)*(-9.0)*0.5)
         if((resn(i).eq.'GLU').and.(atnam(i)(2:4).eq.'OE2'))
     1   sfe(i)=(aburied(i)*(-37.)*0.5)+(aburied(i)*(-9.0)*0.5)
        endif
c       tyrosine
        if((pHvalu.ge.10.1).and.(sasafold(i).gt.thresh))then
         if(atnam(i)(2:4).eq.'OH ')sfe(i)=aburied(i)*(-37.)
        endif
c       cterminal carboxylate
        if((pHvalu.ge.3.0).and.(sasafold(i).gt.thresh))then
         if(atnam(i)(2:4).eq.'OXT')sfe(i)=aburied(i)*(-37.0)
         if(atnam(i)(2:4).eq.'OX ')sfe(i)=aburied(i)*(-37.0)
         if(atnam(i)(2:4).eq.'OT ')sfe(i)=aburied(i)*(-37.0)
         if(atnam(i)(2:4).eq.'OT2')sfe(i)=aburied(i)*(-37.0)
        endif
c       N-terminal amino
c       Rely on N-terminal amino acids to be called resid 1.
        if(resid(i).eq.1)then
         if((pHvalu.le.8.0).and.(sasafold(i).gt.thresh))then
          if(atnam(i)(2:4).eq.'N  ')then
c          print*, 'I found charged N-term resid',resid(i)
           sfe(i)=aburied(i)*(-38.0)
          endif
         endif
       endif
       sumnrgy=sumnrgy+sfe(i)
c      dividing the solvation energy by 10 to fit within columns allowed for Bfactor
       write(8,92)atnam(i),resn(i),ch(i),resid(i),x(i),
     1  y(i),z(i),aburied(i),sfe(i)/10.0
92     format('ATOM      1 ',a5,a3,1x,a1,i4,4x,3f8.3,f6.1,f6.1)
91     continue
       write(6,15)chainlist(j),sumnrgy,pHvalu,thresh
15     format('Solvation free energy for chain ',a1,' is',
     1 f9.1,' cal/mol at pH',f5.1,' assuming charged if sasafold.gt.',
     1 f5.1)
93    continue
c
c     kluge KLUGE kluge KLUGE kluge KLUGE kluge KLUGE kluge
c     kluge KLUGE kluge KLUGE kluge KLUGE kluge KLUGE kluge
c     Calculate solvation free energy, using the rules outlined in the remarks section
c     for this particular kluge version.
c
c     Set default charge flag to 1 for all atoms, meaning they are charged (guilty until proven innocent). 
c     charge flag = 0 means the charge is paired.
      do 145, i=1,n
       charge(i)=1
       asp(i)=-999
145   continue
c     Find all atoms with capability of positive charge (His, Lys, Arg)
      do 155, i=1,n
       if(((resn(i).eq.'HIS').and.
     1     ((atnam(i)(2:4).eq.'ND1').or.
     1      (atnam(i)(2:4).eq.'NE2'))).or.
     1    ((resn(i).eq.'LYS').and.
     1      (atnam(i)(2:4).eq.'NZ ')).or.
     1    ((resn(i).eq.'ARG').and.
     1     ((atnam(i)(2:4).eq.'NH1').or.
     1      (atnam(i)(2:4).eq.'NH2').or.
     1      (atnam(i)(2:4).eq.'NE '))).or.
     1    (resn(i).eq.'NTE'))then
c       Find all atoms with capability of negative charge  (Asp, Glu)
        do 154, ii=1,n
        if(((resn(ii).eq.'ASP').and.
     1      ((atnam(ii)(2:4).eq.'OD1').or.
     1       (atnam(ii)(2:4).eq.'OD2'))).or.
     1     ((resn(ii).eq.'GLU').and.
     1      ((atnam(ii)(2:4).eq.'OE1').or.
     1       (atnam(ii)(2:4).eq.'OE2'))).or.
     1      (resn(ii).eq.'CTE'))then
         dionx=(x(i)-x(ii))**2
         diony=(y(i)-y(ii))**2
         dionz=(z(i)-z(ii))**2
         dist=sqrt(dionx+diony+dionz)
c        If distance between atoms with complementary charges < 4.5 Angstroms
c        **AND** sasafold of each atom < 50 A^2, then print
c        and set charge state for both atoms in ion pair from one to zero.
         if(dist.le.4.5)then
          if((sasafold(i).lt.50.00).and.(sasafold(ii).lt.50.00))then
           charge(i)=0
           charge(ii)=0
           disttemp=dist
           if(dist.lt.2.8)disttemp=2.8
           asptemp=-9.0*(((disttemp-2.8)/2.8)**2)
           if(asptemp.gt.asp(i) )asp(i) =asptemp
           if(asptemp.gt.asp(ii))asp(ii)=asptemp
           write(6,152)ch(i),resid(i),resn(i),atnam(i)(2:4),
     1     charge(i),sasafold(i),asp(i)
           write(6,153)ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1     charge(ii),sasafold(ii),asp(ii),dist
152        format( 'PAIRED CHARG+ ',a1,i4,1x,a3,1x,a3,' chrg=',f5.1,
     1     ' sasa=',f5.1,' a.s.p.=',f4.1)
153        format( 'PAIRED CHARG- ',a1,i4,1x,a3,1x,a3,' chrg=',f5.1,
     1     ' sasa=',f5.1,' a.s.p.=',f4.1,' distance=',f5.1)
          endif
         endif
        endif
154    continue
       endif
155   continue
c     End of complementary charge search

C     If charged amino acid is binuclear (i.e. Asp or Glu) then...
c     If one oxygen is in an ion pair, set charge(i) for both oxygens to zero.
c     If neither oxygen is in an ion pair, set charge(i) to zero for the atom that is least buried.
      do 255, i=1,n
c      Glu search for binuclear mates
c      Glu search for binuclear mates
c      Find OE1                       
       if((resn(i).eq.'GLU').and.(atnam(i)(2:4).eq.'OE1'))then
c       Find mated nucleus OE2
        do 254, ii=1,n
         if((resn(ii).eq.'GLU').and.(atnam(ii)(2:4).eq.'OE2').and.
     1    (resid(ii).eq.resid(i)).and.(ch(ii).eq.ch(i)))then
c         If OE1 is in an ion pair, set charge(i) for OE2 to zero.
          if((charge(i).eq.0).and.(charge(ii).eq.1))then
           charge(ii)=0
           write(6,253)ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1     charge(ii),ch(i),resid(i),resn(i),atnam(i)(2:4),charge(i)
253        format( 'BINUC NEUTRAL ',a1,i4,1x,a3,1x,a3,' chrg=',f5.1,
     1     ' because binuclear mate ',a1,i4,1x,a3,1x,a3,' chrg=',f5.1,
     1     ' is in ionpair')
           goto 255
          endif
c         If OE2 is in an ion pair, set charge(i) for OE1 to zero.
          if((charge(i).eq.1).and.(charge(ii).eq.0))then
           charge(i)=0
           write(6,253)ch(i),resid(i),resn(i),atnam(i)(2:4),
     1     charge(i),ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1     charge(ii)
           goto 255
          endif
c         If neither OE1 or OE2 is in an ion pair, set charge(i)=0 for least buried oxygen.
          if((charge(i).eq.1).and.(charge(ii).eq.1))then
c          If OE1 is least buried oxygen, it is neutralized
           if(aburied(i).le.aburied(ii))then
            charge(ii)=0
            write(6,256)ch(i),resid(i),resn(i),atnam(i)(2:4),charge(i),
     1      aburied(i)
            write(6,257)ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1      charge(ii),aburied(ii)
256         format( 'BINUC CHARGE- ',a1,i4,1x,a3,1x,a3,' chrg=',f5.1,
     1      ' abu=',f5.1,' less buried part of binuclear acid unpaired')
257         format( 'BINUC CHARGE0 ',a1,i4,1x,a3,1x,a3,' chrg=',f5.1,
     1      ' abu=',f5.1,' more buried part of binuclear acid unpaired')
            goto 255
           endif
c          If OE2 is least buried oxygen, it is neutralized.
           if(aburied(ii).lt.aburied(i))then
            charge(i)=0
            write(6,256)ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1      charge(ii),aburied(ii)
            write(6,257)ch(i),resid(i),resn(i),atnam(i)(2:4),charge(i),
     1      aburied(i)
            goto 255
           endif
          endif
         endif
254     continue
       endif
c      Asp search for binuclear mates
c      Asp search for binuclear mates
c      Find OD1                       
       if((resn(i).eq.'ASP').and.(atnam(i)(2:4).eq.'OD1'))then
c       Find mated nucleus OD2
        do 264, ii=1,n
         if((resn(ii).eq.'ASP').and.(atnam(ii)(2:4).eq.'OD2').and.
     1    (resid(ii).eq.resid(i)).and.(ch(ii).eq.ch(i)))then
c         If OD1 is in an ion pair, set charge(i) for OD2 to zero.
          if((charge(i).eq.0).and.(charge(ii).eq.1))then
           charge(ii)=0
           write(6,253)ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1     charge(ii),ch(i),resid(i),resn(i),atnam(i)(2:4),charge(i)
           goto 255
          endif
c         If OD2 is in an ion pair, set charge(i) for OD1 to zero.
          if((charge(i).eq.1).and.(charge(ii).eq.0))then
           charge(i)=0
           write(6,253)ch(i),resid(i),resn(i),atnam(i)(2:4),
     1     charge(i),ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1     charge(ii)
           goto 255
          endif
c         If neither OD1 or OD2 is in an ion pair, set charge(i)=0 for least buried oxygen.
          if((charge(i).eq.1).and.(charge(ii).eq.1))then
c          If OD1 is least buried oxygen, it is neutralized
           if(aburied(i).le.aburied(ii))then
            charge(ii)=0
            write(6,256)ch(i),resid(i),resn(i),atnam(i)(2:4),charge(i),
     1      aburied(i)
            write(6,257)ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1      charge(ii),aburied(ii)
            goto 255
           endif
c          If OD2 is least buried oxygen, it is neutralized.
           if(aburied(ii).lt.aburied(i))then
            charge(i)=0
            write(6,256)ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1      charge(ii),aburied(ii)
            write(6,257)ch(i),resid(i),resn(i),atnam(i)(2:4),charge(i),
     1      aburied(i)
            goto 255
           endif
          endif
         endif
264     continue
       endif
c      CTE search for binuclear mates
c      CTE search for binuclear mates
c      Find CTE                       
       if((resn(i).eq.'CTE').and.(atnam(i)(2:4).eq.'O  '))then
c       Find mated nucleus OXT
        do 274, ii=1,n
         if((resn(ii).eq.'CTE').and.(atnam(ii)(2:4).eq.'OXT').and.
     1    (resid(ii).eq.resid(i)).and.(ch(ii).eq.ch(i)))then
c         If O   is in an ion pair, set charge(i) for OXT to zero.
          if((charge(i).eq.0).and.(charge(ii).eq.1))then
           charge(ii)=0
           write(6,253)ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1     charge(ii),ch(i),resid(i),resn(i),atnam(i)(2:4),charge(i)
           goto 255
          endif
c         If OXT is in an ion pair, set charge(i) for O   to zero.
          if((charge(i).eq.1).and.(charge(ii).eq.0))then
           charge(i)=0
           write(6,253)ch(i),resid(i),resn(i),atnam(i)(2:4),
     1     charge(i),ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1     charge(ii)
           goto 255
          endif
c         If neither O   or OXT is in an ion pair, set charge(i)=0 for least buried oxygen.
          if((charge(i).eq.1).and.(charge(ii).eq.1))then
c          If O   is least buried oxygen, it is neutralized
           if(aburied(i).le.aburied(ii))then
            charge(ii)=0
            write(6,256)ch(i),resid(i),resn(i),atnam(i)(2:4),charge(i),
     1      aburied(i)
            write(6,257)ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1      charge(ii),aburied(ii)
            goto 255
           endif
c          If OXT is least buried oxygen, it is neutralized.
           if(aburied(ii).lt.aburied(i))then
            charge(i)=0
            write(6,256)ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1      charge(ii),aburied(ii)
            write(6,257)ch(i),resid(i),resn(i),atnam(i)(2:4),charge(i),
     1      aburied(i)
            goto 255
           endif
          endif
         endif
274     continue
       endif
c      His search for binuclear mates
c      His search for binuclear mates
c      Find ND1                       
       if((resn(i).eq.'HIS').and.(atnam(i)(2:4).eq.'ND1'))then
c       Find mated nucleus NE2
        do 284, ii=1,n
         if((resn(ii).eq.'HIS').and.(atnam(ii)(2:4).eq.'NE2').and.
     1    (resid(ii).eq.resid(i)).and.(ch(ii).eq.ch(i)))then
c         If ND1 is in an ion pair, set charge(i) for NE2 to zero.
          if((charge(i).eq.0).and.(charge(ii).eq.1))then
           charge(ii)=0
           write(6,253)ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1     charge(ii),ch(i),resid(i),resn(i),atnam(i)(2:4),charge(i)
           goto 255
          endif
c         If NE2 is in an ion pair, set charge(i) for ND1 to zero.
          if((charge(i).eq.1).and.(charge(ii).eq.0))then
           charge(i)=0
           write(6,253)ch(i),resid(i),resn(i),atnam(i)(2:4),
     1     charge(i),ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1     charge(ii)
           goto 255
          endif
c         If neither ND1 or NE2 is in an ion pair, set charge(i)=0 for least buried nitrogen.
          if((charge(i).eq.1).and.(charge(ii).eq.1))then
c          If ND1 is less buried nitrogen, it is neutralized
           if(aburied(i).le.aburied(ii))then
            charge(ii)=0
            write(6,258)ch(i),resid(i),resn(i),atnam(i)(2:4),charge(i),
     1      aburied(i)
            write(6,259)ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1      charge(ii),aburied(ii)
258         format( 'BINUC CHARGE+ ',a1,i4,1x,a3,1x,a3,' chrg=',f5.1,
     1      ' abu=',f5.1,' less buried part of binuclear base unpaired')
259         format( 'BINUC CHARGE0 ',a1,i4,1x,a3,1x,a3,' chrg=',f5.1,
     1      ' abu=',f5.1,' more buried part of binuclear base unpaired')
            goto 255
           endif
c          If NE2 is less buried nitrogen, it is neutralized.
           if(aburied(ii).lt.aburied(i))then
            charge(i)=0
            write(6,258)ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1      charge(ii),aburied(ii)
            write(6,259)ch(i),resid(i),resn(i),atnam(i)(2:4),charge(i),
     1      aburied(i)
            goto 255
           endif
          endif
         endif
284     continue
       endif
c      Arg search for trinuclear mates
c      Arg search for trinuclear mates
c      Find NE                        
       if((resn(i).eq.'ARG').and.(atnam(i)(2:4).eq.'NE '))then
c       Find mated nucleus NH1
        do 294, ii=1,n
         if((resn(ii).eq.'ARG').and.(atnam(ii)(2:4).eq.'NH1').and.
     1    (resid(ii).eq.resid(i)).and.(ch(ii).eq.ch(i)))goto 293
294     continue
c       Find mated nucleus NH2
293     do 295, iii=1,n
         if((resn(iii).eq.'ARG').and.(atnam(iii)(2:4).eq.'NH2').and.
     1    (resid(iii).eq.resid(i)).and.(ch(iii).eq.ch(i)))then
c         If NE  is in an ion pair, set charge for NH1 and NH2 to zero.
          if(charge(i).eq.0)then
           charge(ii)=0
           charge(iii)=0
           write(6,260)ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1     charge(ii),ch(i),resid(i),resn(i),atnam(i)(2:4),charge(i)
           write(6,260)ch(iii),resid(iii),resn(iii),atnam(iii)(2:4),
     1     charge(iii),ch(i),resid(i),resn(i),atnam(i)(2:4),charge(i)
260        format( 'TRNUC NEUTRAL ',a1,i4,1x,a3,1x,a3,' chrg=',f5.1,
     1     ' bcause trinuclear mate ',a1,i4,1x,a3,1x,a3,' chrg=',f5.1,
     1     ' is in ionpair')
           goto 255
          endif
c         If NH1 is in an ion pair, set charge for NE and NH2 to zero.
          if(charge(ii).eq.0)then
           charge(i)=0
           charge(iii)=0
           write(6,260)ch(i),resid(i),resn(i),atnam(i)(2:4),
     1     charge(i),ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),charge(ii)
           write(6,260)ch(iii),resid(iii),resn(iii),atnam(iii)(2:4),
     1     charge(iii),ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1     charge(ii)
           goto 255
          endif
c         If NH2 is in an ion pair, set charge for NE and NH1 to zero.
          if(charge(iii).eq.0)then
           charge(i)=0
           charge(ii)=0
           write(6,260)ch(i),resid(i),resn(i),atnam(i)(2:4),
     1     charge(i),ch(iii),resid(iii),resn(iii),atnam(iii)(2:4),
     1     charge(iii)
           write(6,260)ch(ii),resid(ii),resn(iii),atnam(ii)(2:4),
     1     charge(ii),ch(iii),resid(iii),resn(iii),atnam(iii)(2:4),
     1     charge(iii)
           goto 255
          endif
c         If neither NE nor NH1 nor NH2 is in an ion pair, set charge=0 for least buried nitrogen.
          if((charge(i).eq.1).and.(charge(ii).eq.1).and.
     1       (charge(iii).eq.1))then
c          If NE is least buried nitrogen, it is neutralized
           if((aburied(i).le.aburied(ii)).and.
     1        (aburied(i).le.aburied(iii)))then
            charge(ii)=0
            charge(iii)=0
            write(6,261)ch(i),resid(i),resn(i),atnam(i)(2:4),
     1      charge(i),aburied(i)
            write(6,262)ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1      charge(ii),aburied(ii)
            write(6,262)ch(iii),resid(iii),resn(iii),atnam(iii)(2:4),
     1      charge(iii),aburied(iii)
261         format( 'TRNUC CHARGE+ ',a1,i4,1x,a3,1x,a3,' chrg=',f5.1,
     1      ' abu=',f5.1,' least buried nitrogen of unpaired Arg')
262         format( 'TRNUC CHARGE0 ',a1,i4,1x,a3,1x,a3,' chrg=',f5.1,
     1      ' abu=',f5.1,' more  buried nitrogen of unpaired Arg')
            goto 255
           endif
c          If NH1 is least buried nitrogen, it is neutralized.
           if((aburied(ii).le.aburied(i)).and.
     1        (aburied(ii).le.aburied(iii)))then
            charge(i)=0
            charge(iii)=0
            write(6,261)ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1      charge(ii),aburied(ii)
            write(6,262)ch(i),resid(i),resn(i),atnam(i)(2:4),
     1      charge(i),aburied(i)
            write(6,262)ch(iii),resid(iii),resn(iii),atnam(iii)(2:4),
     1      charge(iii),aburied(iii)
            goto 255
           endif
c          If NH2 is least buried nitrogen, it is neutralized.
           if((aburied(iii).le.aburied(i)).and.
     1        (aburied(iii).le.aburied(ii)))then
            charge(i)=0
            charge(ii)=0
            write(6,261)ch(iii),resid(iii),resn(iii),atnam(iii)(2:4),
     1      charge(iii),aburied(iii)
            write(6,262)ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),
     1      charge(ii),aburied(ii)
            write(6,262)ch(i),resid(i),resn(i),atnam(i)(2:4),
     1      charge(i),aburied(i)
            goto 255
           endif
          endif
         endif
295     continue
       endif
255   continue

c
c     Begin search for H-bonds
c     Set default Hbond count to 0 for all atoms.
      do 146, i=1,n
       nhbond(i)=0
146   continue
c     Write the load pdbfile record to pymol file
      write (10,160)pdbinfile1
160   format('load ',a90,', proteorella')
c     Find all potential H-bond acceptors, obligate or otherwise
      do 156, i=1,n
       if(atnam(i)(2:2).eq.'O')then
c       Find all potential H-bond donors, obligate or otherwise
        do 157, ii=1,n
        if((atnam(ii)(2:2).eq.'N').or.
     1     (atnam(ii)(2:4).eq.'OG ').or.
     1     (atnam(ii)(2:4).eq.'OH ').or.
     1     (atnam(ii)(2:4).eq.'OG1'))then
         dionx=(x(i)-x(ii))**2
         diony=(y(i)-y(ii))**2
         dionz=(z(i)-z(ii))**2
         dist=sqrt(dionx+diony+dionz)
c        If distance between atoms with complementary charges < 3.2 A and >2.3 A, then 
c        write the distance to a pymol file.
         if((dist.lt.3.3).and.(dist.ge.2.3))then
          if((ch(i).eq.ch(ii)).and.(resid(i).eq.resid(ii)))goto 157
          write(10,158)ch(i),resid(i),
     1    atnam(i)(2:4),ch(ii),resid(ii),atnam(ii)(2:4)
158       format('distance (chain ',a1,' and resid ',i4,' and name ',
     1    a3,'), (chain ',a1,' and resid ',i4,' and name ',a3,')')
c         write(6,159)ch(i),resid(i),resn(i),
c    1    atnam(i)(2:4),ch(ii),resid(ii),resn(ii),atnam(ii)(2:4),dist
159       format( 'HBOND FOUND   ',a1,i4,1x,a3,1x,a3,5x,
     1    a1,i4,1x,a3,1x,a3,' distance=',f5.1)
          nhbond(i)=nhbond(i)+1
          nhbond(ii)=nhbond(ii)+1
         endif
        endif
157    continue
       endif
156   continue
c     End of hbond search
c
c     We insert here the routines needed to evaluate entropic contributions.
c     We want to scale the entropy of each amino acid type by the percentage of side chain surface area 
c     that is buried in the assembly. Begin by evaluating the surface area buried per residue.
c     Index area values by residue count, rescnt.
      rescnt=0
      do 401,i=1,n
c      Skip main chain atoms because these entropy terms pertain only to side chains.
       if(atnam(i).eq.' N  ')goto 401
       if(atnam(i).eq.' C  ')goto 401
       if(atnam(i).eq.' O  ')goto 401
       if((atnam(i).eq.' CA ').and.(resn(i).ne.'GLY'))goto 401
c      If this atom is glycine, substitute the side chain area with alpha carbon area.
       if((atnam(i).eq.' CA ').and.(resn(i).eq.'GLY'))then
        rescnt=rescnt+1
        nthresresid(rescnt)=resid(i)
        nthresch(rescnt)=ch(i)
        nthressasafold(rescnt)=sasafold(i)
        nthressasaisol(rescnt)=sasaisol(i)
        goto 401
       endif
c      If this atom is the first atom encountered, record SASAs.
       if(rescnt.eq.0)then
        rescnt=1
        nthresresid(rescnt)=resid(i)
        nthresch(rescnt)=ch(i)
        nthressasafold(rescnt)=sasafold(i)
        nthressasaisol(rescnt)=sasaisol(i)
        goto 401
       endif
       do 402, m=1,rescnt 
c       If this atom matches the chain ID and resid of a previous atom, 
c       record SASAs for this atom as a function of rescnt. Proceed to next atom.
        if((resid(i).eq.nthresresid(m)).and.(ch(i).eq.nthresch(m)))then
         nthressasafold(rescnt)=nthressasafold(rescnt)+sasafold(i)
         nthressasaisol(rescnt)=nthressasaisol(rescnt)+sasaisol(i)
         goto 401
        endif
402    continue
c      If this atom does not match previous atom chain ID and resid, 
c      record SASAs for this atom as a function of rescnt.
       rescnt=rescnt+1
       nthresresid(rescnt)=resid(i)
       nthresch(rescnt)=ch(i)
       nthressasafold(rescnt)=sasafold(i)
       nthressasaisol(rescnt)=sasaisol(i)
401   continue
      print*,rescnt,' residues total among all chains.'
c
c     Calculate the percent area buried per residue.
      do 403, m=1,rescnt
       tempabu=nthressasaisol(m)-nthressasafold(m)
       if(tempabu.lt.0)temp=0
       nthrespctabu(m)=tempabu/nthressasaisol(m)
c      print*,m,nthressasaisol(m),nthressasafold(m),nthrespctabu(m)
403   continue
c     Assign the percent area buried per residue to the atoms it comprises.
c     Need to do this to report percent area in B-factor column.
      do 405, i=1,n
       do 406, m=1,rescnt
        if((resid(i).eq.nthresresid(m)).and.(ch(i).eq.nthresch(m)))then
         nthatompctabu(i)=nthrespctabu(m)
         goto 405
        endif
406    continue
405   continue
c     Report percent area buried per residue in pctsideburied.pdb.
      do 407, i=1,n
       write(12,49)i,atnam(i),resn(i),ch(i),resid(i),
     1 x(i),y(i),z(i),nthatompctabu(i)
c      Unsatisfied HB donor and acceptors written to hbunsatisfied.pdb
       if((ch(i).eq.'A').and.(nhbond(i).eq.0).and.
     1 ((atnam(i)(2:2).eq.'N').or.(atnam(i)(2:2).eq.'O')))then
        write(14,52)i,atnam(i),resn(i),ch(i),resid(i),
     1  x(i),y(i),z(i),float(nhbond(i)),sasafold(i)
52      format('ATOM',i7,1x,a5,a3,1x,a1,i4,4x,3f8.3,2f6.2)
       endif
407   continue
c
c     Count the number of each amino acid's occurrence in each chain.
      do 315, j=1,nchains
       nuA(j)=0
       nuC(j)=0
       nuD(j)=0
       nuE(j)=0
       nuF(j)=0
       nuG(j)=0
       nuH(j)=0
       nuI(j)=0
       nuK(j)=0
       nuL(j)=0
       nuM(j)=0
       nuN(j)=0
       nuP(j)=0
       nuQ(j)=0
       nuR(j)=0
       nuS(j)=0
       nuT(j)=0
       nuV(j)=0
       nuW(j)=0
       nuY(j)=0
       do 316, i=1,n
c       Regard one chain at a time.
        if(ch(i).ne.chainlist(j))goto 316
c       TdS is Koehl and Delarue 1994.
c       Scale TdS by the percentage of the side chain surface area that is buried.
c       scaleTdS=1.2*(nthatompctabu(i)*nthatompctabu(i))
        scaleTdS=1.0*(nthatompctabu(i))
        if(resn(i).eq.'ALA')entpatom(i)=scaleTdS*( 0.00)/5.0
        if(resn(i).eq.'CYS')entpatom(i)=scaleTdS*(-0.63)/6.0
        if(resn(i).eq.'ASP')entpatom(i)=scaleTdS*(-0.65)/8.0
        if(resn(i).eq.'GLU')entpatom(i)=scaleTdS*(-1.31)/9.0
        if(resn(i).eq.'PHE')entpatom(i)=scaleTdS*(-0.65)/11.0
        if(resn(i).eq.'GLY')entpatom(i)=scaleTdS*( 0.00)/4.0
        if(resn(i).eq.'HIS')entpatom(i)=scaleTdS*(-0.92)/10.0
        if(resn(i).eq.'ILE')entpatom(i)=scaleTdS*(-0.94)/8.0
        if(resn(i).eq.'LYS')entpatom(i)=scaleTdS*(-1.63)/9.0
        if(resn(i).eq.'LEU')entpatom(i)=scaleTdS*(-0.94)/8.0
        if(resn(i).eq.'MET')entpatom(i)=scaleTdS*(-1.24)/8.0
        if(resn(i).eq.'ASN')entpatom(i)=scaleTdS*(-0.75)/8.0
        if(resn(i).eq.'PRO')entpatom(i)=scaleTdS*(-0.30)/7.0
        if(resn(i).eq.'GLN')entpatom(i)=scaleTdS*(-1.29)/9.0
        if(resn(i).eq.'ARG')entpatom(i)=scaleTdS*(-1.21)/11.0
        if(resn(i).eq.'SER')entpatom(i)=scaleTdS*(-0.43)/6.0
        if(resn(i).eq.'THR')entpatom(i)=scaleTdS*(-0.57)/7.0
        if(resn(i).eq.'VAL')entpatom(i)=scaleTdS*(-0.62)/7.0
        if(resn(i).eq.'TRP')entpatom(i)=scaleTdS*(-1.14)/14.0
        if(resn(i).eq.'TYR')entpatom(i)=scaleTdS*(-1.07)/12.0
c       print*, i,ch(i),sasafold(i),scaleTdS,entpatom(i)
        if(atnam(i).ne.' CA  ')goto 316
        if(resn(i).eq.'ALA')nuA(j)=nuA(j)+1
        if(resn(i).eq.'CYS')nuC(j)=nuC(j)+1
        if(resn(i).eq.'ASP')nuD(j)=nuD(j)+1
        if(resn(i).eq.'GLU')nuE(j)=nuE(j)+1
        if(resn(i).eq.'PHE')nuF(j)=nuF(j)+1
        if(resn(i).eq.'GLY')nuG(j)=nuG(j)+1
        if(resn(i).eq.'HIS')nuH(j)=nuH(j)+1
        if(resn(i).eq.'ILE')nuI(j)=nuI(j)+1
        if(resn(i).eq.'LYS')nuK(j)=nuK(j)+1
        if(resn(i).eq.'LEU')nuL(j)=nuL(j)+1
        if(resn(i).eq.'MET')nuM(j)=nuM(j)+1
        if(resn(i).eq.'ASN')nuN(j)=nuN(j)+1
        if(resn(i).eq.'PRO')nuP(j)=nuP(j)+1
        if(resn(i).eq.'GLN')nuQ(j)=nuQ(j)+1
        if(resn(i).eq.'ARG')nuR(j)=nuR(j)+1
        if(resn(i).eq.'SER')nuS(j)=nuS(j)+1
        if(resn(i).eq.'THR')nuT(j)=nuT(j)+1
        if(resn(i).eq.'VAL')nuV(j)=nuV(j)+1
        if(resn(i).eq.'TRP')nuW(j)=nuW(j)+1
        if(resn(i).eq.'TYR')nuY(j)=nuY(j)+1
316    continue
       write(6,317)nuA(j),nuC(j),nuD(j),nuE(j),nuF(j),nuG(j),nuH(j),
     1      nuI(j),nuK(j),nuL(j),nuM(j),nuN(j),nuP(j),nuQ(j),nuR(j),
     1      nuS(j),nuT(j),nuV(j),nuW(j),nuY(j)
317    format(i3,'A ',i3'C ',i3,'D ',i3,'E ',i3,'F ',i3,'G ',i3'H ',
     1        i3,'I ',i3'K ',i3,'L ',i3,'M ',i3,'N ',i3,'P ',i3'Q ',
     1        i3,'R ',i3'S ',i3,'T ',i3,'V ',i3,'W ',i3,'Y')
315   continue
c     Calculate the entropic penalty of each chain.
c     TdS is mean of 6 models, Doig and Sternberg 1995 Protein Sci. 4, 2247.
      do 318, j=1,nchains
       entch(j)=0
       do 320, i=1,n
c       Regard one chain at a time.
        if(ch(i).ne.chainlist(j))goto 320
        entch(j)=entch(j)+entpatom(i)
320    continue
       write(6,319)chainlist(j),-entch(j)
319    format('Entropic penalty T*dS for chain',a2,' =',f8.2)
318   continue
c
c     Final big do loop to calculate SFE from the atomic solvation parameters formulated so carefully above.
c     And combine with entropy terms and write the results to pdb file.
      do 193, j=1,nchains
       sumnrgy=0.0
       nHB__satisf=0
       nHBunsatisf=0
       nHBunsatisfbur=0
       do 191, i=1,n
c       Regard one chain at a time.
        if(ch(i).ne.chainlist(j))goto 191
        aburied(i)=sasaisol(i)-sasafold(i)
c       CARBON
c       CARBON
        if(atnam(i)(2:2).eq.'C')then
         asp(i)=18.0
         sfe(i)=aburied(i)*asp(i)
        endif
c       NITROGEN
c       NITROGEN
        if(atnam(i)(2:2).eq.'N')then
         if(nhbond(i).eq.0)then
          nHBunsatisf=nHBunsatisf+1
          if(sasafold(i).le.10)nHBunsatisfbur=nHBunsatisfbur+1
         endif
         if(nhbond(i).gt.0)nHB__satisf=nHB__satisf+1
c        Set uncharged A.S.P. as default value.
         if(asp(i).eq.-999)asp(i)=-9.0
         sfe(i)=aburied(i)*asp(i)
c        if N atom is in the backbone and is hydrogen bonded, set A.S.P. to zero
         if((atnam(i)(2:4).eq.'N  ').and.(nhbond(i).gt.0))then
          asp(i)=0.0
          sfe(i)=aburied(i)*asp(i)
         endif
c        if N atom is in the N-terminus
         if(resn(i).eq.'NTE')then
          if(charge(i).eq.0)asp(i)=asp(i)
          if(charge(i).eq.1)asp(i)=-38
          if(pHvalu.gt.8)asp(i)=-9
          sfe(i)=aburied(i)*asp(i)
         endif
c        if N atom is in Lys sidechain 
         if((atnam(i)(2:4).eq.'NZ ').and.(resn(i).eq.'LYS'))then
          if(charge(i).eq.0)asp(i)=asp(i)
          if(charge(i).eq.1)asp(i)=-38
          if(pHvalu.gt.10)asp(i)=-9
          sfe(i)=aburied(i)*asp(i)
         endif
c        if N atom is in Arg sidechain 
         if(((atnam(i)(2:4).eq.'NE ').or.
     1       (atnam(i)(2:4).eq.'NH1').or.
     1       (atnam(i)(2:4).eq.'NH2')).and.
     1       (resn(i).eq.'ARG'))then
          if(charge(i).eq.0)asp(i)=asp(i)
          if(charge(i).eq.1)asp(i)=-38
          if(pHvalu.gt.10.5)asp(i)=-9
          sfe(i)=aburied(i)*asp(i)
         endif
c        if N atom is in His sidechain 
         if(((atnam(i)(2:4).eq.'ND1').or.
     1       (atnam(i)(2:4).eq.'NE2')).and.
     1       (resn(i).eq.'HIS'))then
          if(charge(i).eq.0)asp(i)=asp(i)
          if(charge(i).eq.1)asp(i)=-38
          if(pHvalu.gt.6.08)asp(i)=-9
          sfe(i)=aburied(i)*asp(i)
         endif
c        if N atom is in Gln sidechain and has 2 hydrogen bonds
c        **AND** sasafold < 5 A^2, then set A.S.P. to zero
         if((atnam(i)(2:4).eq.'NE2').and.(nhbond(i).gt.0).and.
     1   (resn(i).eq.'GLN').and.(sasafold(i).lt.5))then
c         Find how many H-bonds are on OE1
          do 301, ii=1,n
           if((resn(ii).eq.'GLN').and.(atnam(ii)(2:4).eq.'OE1').and.
     1     (resid(ii).eq.resid(i)).and.(ch(ii).eq.ch(i)))then
            if(nhbond(ii).gt.0)then
             if(sasafold(ii).lt.5)then
              asp(i)=0.0
              write(6,331)resn(i),ch(i),resid(i),atnam(i)(2:4),
     1        nhbond(i),sasafold(i),resn(ii),ch(ii),resid(ii),
     1        atnam(ii)(2:4),nhbond(ii),sasafold(ii)
331           format('Buried ',a3,1x,a1,1x,i4,1x,a3,
     1        ' has',i2,' hbonds, sasa=',f6.1,
     1        ' conspiring with ',a3,1x,a1,1x,i4,1x,a3,
     1        ' has',i2,' hbonds, sasa=',f6.1)
             endif
            endif
            goto 302
           endif
301       continue
302       sfe(i)=aburied(i)*asp(i)
         endif
c        if N atom is in Asn sidechain and has 2 hydrogen bonds 
c        **AND** sasafold < 5 A^2, then set A.S.P. to zero
         if((atnam(i)(2:4).eq.'ND2').and.(nhbond(i).gt.0).and.
     1   (resn(i).eq.'ASN').and.(sasafold(i).lt.5))then
c         Find how many H-bonds are on OD1
          do 303, ii=1,n
           if((resn(ii).eq.'ASN').and.(atnam(ii)(2:4).eq.'OD1').and.
     1     (resid(ii).eq.resid(i)).and.(ch(ii).eq.ch(i)))then
            if(nhbond(ii).gt.0)then
             if(sasafold(ii).lt.5)then
              asp(i)=0.0
              write(6,331)resn(i),ch(i),resid(i),atnam(i)(2:4),
     1        nhbond(i),sasafold(i),resn(ii),ch(ii),resid(ii),
     1        atnam(ii)(2:4),nhbond(ii),sasafold(ii)
             endif
            endif
            goto 304
           endif
303       continue
304       sfe(i)=aburied(i)*asp(i)
         endif
c        if NXT atom is amidated Cterm (ATE) and has 2 hydrogen bonds 
c        then set A.S.P. to zero. 
c        No requirement for sasafold < 5 A^2 because these are backbone atoms
         if((atnam(i)(2:4).eq.'NXT').and.(nhbond(i).gt.0).and.
     1   (resn(i).eq.'ATE'))then
c         Find how many H-bonds are on O
          do 323, ii=1,n
           if((atnam(ii)(2:4).eq.'O  ').and.
     1     (resid(ii).eq.resid(i)).and.(ch(ii).eq.ch(i)))then
            if(nhbond(ii).gt.0)then
              asp(i)=0.0
              write(6,331)resn(i),ch(i),resid(i),atnam(i)(2:4),
     1        nhbond(i),sasafold(i),resn(ii),ch(ii),resid(ii),
     1        atnam(ii)(2:4),nhbond(ii),sasafold(ii)
            endif
            goto 324
           endif
323       continue
324       sfe(i)=aburied(i)*asp(i)
         endif

        endif
c       OXYGEN  
c       OXYGEN  
        if(atnam(i)(2:2).eq.'O')then
         if(nhbond(i).eq.0)then
          nHBunsatisf=nHBunsatisf+1
          if(sasafold(i).le.10)nHBunsatisfbur=nHBunsatisfbur+1
         endif
         if(nhbond(i).gt.0)nHB__satisf=nHB__satisf+1
c        Set uncharged A.S.P. as default value.
         if(asp(i).eq.-999)asp(i)=-9.0
         sfe(i)=aburied(i)*asp(i)
c        if O atom is in the backbone and is hydrogen bonded, set A.S.P. to zero
         if((atnam(i)(2:4).eq.'O  ').and.(nhbond(i).gt.0))then
          asp(i)=0.0
          sfe(i)=aburied(i)*asp(i)
         endif
c        if O atom is in the C-terminus
         if(resn(i).eq.'CTE')then
          if(charge(i).eq.0)asp(i)=asp(i)
          if(charge(i).eq.1)asp(i)=-37
          if(pHvalu.lt.3.0)asp(i)=-9
          sfe(i)=aburied(i)*asp(i)
         endif
c        if O atom is in Tyr sidechain 
         if((atnam(i)(2:4).eq.'OH ').and.(resn(i).eq.'TYR'))then
          if(charge(i).eq.0)asp(i)=asp(i)
          if(charge(i).eq.1)asp(i)=-37
          if(pHvalu.lt.10.1)asp(i)=-9
          sfe(i)=aburied(i)*asp(i)
         endif
c        if O atom is in Asp sidechain 
         if(((atnam(i)(2:4).eq.'OD1').or.
     1       (atnam(i)(2:4).eq.'OD2')).and.
     1       (resn(i).eq.'ASP'))then
          if(charge(i).eq.0)asp(i)=asp(i)
          if(charge(i).eq.1)asp(i)=-37
          if(pHvalu.lt.3.9)asp(i)=-9
          sfe(i)=aburied(i)*asp(i)
         endif
c        if O atom is in Glu sidechain 
         if(((atnam(i)(2:4).eq.'OE1').or.
     1       (atnam(i)(2:4).eq.'OE2')).and.
     1       (resn(i).eq.'GLU'))then
          if(charge(i).eq.0)asp(i)=asp(i)
          if(charge(i).eq.1)asp(i)=-37
          if(pHvalu.lt.4.3)asp(i)=-9
          sfe(i)=aburied(i)*asp(i)
         endif
c        if O atom is in Gln sidechain and has 2 hydrogen bonds 
c        **AND** sasafold < 5 A^2, then set A.S.P. to zero
         if((atnam(i)(2:4).eq.'OE1').and.(nhbond(i).gt.0).and.
     1   (resn(i).eq.'GLN').and.(sasafold(i).lt.5))then
c         Find how many H-bonds are on NE2
          do 305, ii=1,n
           if((resn(ii).eq.'GLN').and.(atnam(ii)(2:4).eq.'NE2').and.
     1     (resid(ii).eq.resid(i)).and.(ch(ii).eq.ch(i)))then
            if(nhbond(ii).gt.0)then
             if(sasafold(ii).lt.5)then
              asp(i)=0.0
              write(6,331)resn(i),ch(i),resid(i),atnam(i)(2:4),
     1        nhbond(i),sasafold(i),resn(ii),ch(ii),resid(ii),
     1        atnam(ii)(2:4),nhbond(ii),sasafold(ii)
             endif
            endif
            goto 306
           endif
305       continue
306       sfe(i)=aburied(i)*asp(i)
         endif
c        if N atom is in Asn sidechain and has 2 hydrogen bonds 
c        **AND** sasafold < 5 A^2, then set A.S.P. to zero
         if((atnam(i)(2:4).eq.'OD1').and.(nhbond(i).gt.0).and.
     1   (resn(i).eq.'ASN').and.(sasafold(i).lt.5))then
c         Find how many H-bonds are on ND2
          do 307, ii=1,n
           if((resn(ii).eq.'ASN').and.(atnam(ii)(2:4).eq.'ND2').and.
     1     (resid(ii).eq.resid(i)).and.(ch(ii).eq.ch(i)))then
            if(nhbond(ii).gt.0)then
             if(sasafold(ii).lt.5)then
              asp(i)=0.0
              write(6,331)resn(i),ch(i),resid(i),atnam(i)(2:4),
     1        nhbond(i),sasafold(i),resn(ii),ch(ii),resid(ii),
     1        atnam(ii)(2:4),nhbond(ii),sasafold(ii)
             endif
            endif
            goto 308
           endif
307       continue
308       sfe(i)=aburied(i)*asp(i)
         endif
        endif
c       SULFUR  
c       SULFUR  
        if(atnam(i)(2:2).eq.'S')then
         asp(i)=-5.0
         sfe(i)=aburied(i)*asp(i)
        endif

       sumnrgy=sumnrgy+sfe(i)
c      Write sfeklugeX.Y.pdb
c      Dividing the solvation energy by 1000 to fit within columns allowed for Bfactor (kcal)
       dG=(sfe(i)/1000.0)+entpatom(i)
       write(9,192)atnam(i),resn(i),ch(i),resid(i),x(i),
     1  y(i),z(i),aburied(i),dG
192     format('ATOM      1 ',a5,a3,1x,a1,i4,4x,3f8.3,f6.2,f6.2)
c      Write aspklugeX.Y.pdb
c      Put A.S.P. in the Bfactor column 
       write(11,192)atnam(i),resn(i),ch(i),resid(i),x(i),
     1  y(i),z(i),aburied(i),asp(i)
191     continue
       sumnrgy=sumnrgy/1000.0
       sumnrgyperresid=sumnrgy/nresid(j)
       sumnrgyandentropy=sumnrgy+entch(j)
       HBunsatperresid=float(nHBunsatisf)/nresid(j)
       HBunsatburperresid=float(nHBunsatisfbur)/nresid(j)
       write(6,115)chainlist(j),nresid(j),-sumnrgy,-sumnrgyperresid,
     1 pHvalu,-entch(j),-sumnrgyandentropy,-sumnrgyandentropy/nresid(j),
     1 nHB__satisf,nHBunsatisf,HBunsatperresid,nHBunsatisfbur,
     1 HBunsatburperresid
115     format('Solvation free energy for chain ',a1,' with nresid= ',
     1 i4,' is ',f6.1,' kcal/mol or', f7.3, ' kcal/aa using kluge 7.2 '
     1 'at pH',f5.1,' TdS=',f8.2,' dG=',f8.2,' dG/res=',f8.2,
     1 ' nHBsatisf=',i5, ' nHBunsatisf=',i5,' HBunsat/residue=',f5.2,
     1 ' nHBunsatburied=',i5,' HBunsatburied/residue=',f5.2)
193   continue
c     End of Final Big do loop. Hurray!
c
c
c      Sum energy per residue then write sum to pdb file, B factor column
c      to color atoms by residue energy. sfeklugeX.Y.sumperresid.pdb
c      If this atom is the first atom encountered, record energy, dG.
      rescnt=0
      do 501, i=1,n
       dG=(sfe(i)/1000.0)+entpatom(i)
       if(rescnt.eq.0)then
        rescnt=1
        nthresresid(rescnt)=resid(i)
        nthresch(rescnt)=ch(i)
        nthresdG(rescnt)=dG
        goto 501
       endif
       do 502, m=1,rescnt
c       If this atom matches the chain ID and resid of a previous atom, 
c       record dG for this atom as a function of rescnt. Proceed to next atom.
        if((resid(i).eq.nthresresid(m)).and.(ch(i).eq.nthresch(m)))then
         nthresdG(rescnt)=nthresdG(rescnt)+dG
         goto 501
        endif
502    continue
c      If this atom does not match previous atom chain ID and resid, 
c      record dG for this atom as a function of rescnt.
       rescnt=rescnt+1
       nthresresid(rescnt)=resid(i)
       nthresch(rescnt)=ch(i)
       nthresdG(rescnt)=dG
c      print*,i,rescnt,nthresdG(rescnt) 
501   continue
c
c
c     Assign the dG per residue to the atoms it comprises.
c     Need to do this to report percent area in B-factor column 
c     of sfeklugeX.Y.sumperresid.pdb
      do 505, i=1,n
       do 506, m=1,rescnt
        if((resid(i).eq.nthresresid(m)).and.(ch(i).eq.nthresch(m)))then
         nthatomdG(i)=nthresdG(m)
         goto 505
        endif
506    continue
505   continue
c     Report percent area buried per residue in sfeklugeX.Y.sumperresid.pdb
      do 507, i=1,n
       write(13,49)i,atnam(i),resn(i),ch(i),resid(i),
     1 x(i),y(i),z(i),nthatomdG(i)
507   continue
      go to 999
c
998   print*,'Missing atoms from ',ch(n-1),resid(n-1),prevresn,
     1 natomprevresid
c
 999	stop
	end
