c        1         2         3         4         5         6         7
c23456789012345678901234567890123456789012345678901234567890123456789012
c***********************************************************************
      program nsolmain
c
c     NSOL, version 2.4
c
c     Numerical calculation program for Molecular Surface, Volume, and
c     Solvation Energy
c
c     Copyright (c) 2002-2003 Masato Masuya
c     All rights reserved.
c
c     Redistribution and use in source and binary forms, with or without
c     modification, are permitted provided that the following conditions
c     are met:
c
c     1. Redistributions of source code must retain the above copyright
c     notice, this list of conditions and the following disclaimer.
c
c     2. Redistributions in binary form must reproduce the above
c     copyright notice, this list of conditions and the following
c     disclaimer in the documentation and/or other materials provided
c     with the distribution.
c
c     3. The name of the author may not be used to endorse or promote
c     products derived from this software without specific prior written
c     permission.
c
c     THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS
c     OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
c     WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
c     ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
c     DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
c     DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
c     GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
c     INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
c     WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
c     NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
c     SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
c
c     CHANGES:
c
c     2003/12/21 v.2.4 SPEED UP
c     2002/03/20 v.2.0 INITIAL VERSION
c
c     This program is optimized for vector processor. If you want to use
c     NSOL on a scalar processor, you may prefer NSOL 1.x.
c
      parameter(maxatom=10000)
      parameter(maxsurf=10000)
      character pdbfile*512, surffile*512, arg*256
      integer mode, natom, i, nsurf, set
      integer index(maxatom,3), rank(maxatom,3)
      integer nls(maxatom,3), nll(maxatom,3), eset
      character atomch(maxatom)*4, resch(maxatom)*3
      double precision coord(maxatom,3), rmax, radi(maxatom)
      double precision surf(maxsurf,3), vradi(maxatom)
      double precision gpara(maxatom), hpara(maxatom), cpara(maxatom)
      real time1, time2, tar1(2), tar2(2)
      double precision area, volume, esol, varea, vvolume, temp
      logical debug
c     
      write(*,'(A)') 
     &     '--< NSOL(F), version 2.4 '
     &     // 'by Masato Masuya(CCC, Kagoshima Univ.) >--'
c
      temp = 298.0d0
      eset = 5
c
      if (iargc().eq.2) then
         call getarg(1,pdbfile)
         call getarg(2,surffile)
      else if (iargc().eq.3) then
         call getarg(1,pdbfile)
         call getarg(2,surffile)
         call getarg(3,arg)
         read(arg,*) temp 
         eset = 1
      else
         write(*,'(A)') 
     &        'Usage: nsol pdbfile surffile [temp]'
         stop
      end if
c
      mode = 1
c     
c     mode = 1 - all atoms
c     mode = 2 - all atoms except hydrogen atoms
c     mode = 3 - main chain only
c
      if (eset.eq.1) then
         write(*,'(a)') 'temp dep.: on'
      else 
         write(*,'(a)') 'temp dep.: off'
      end if
c
      call surfread(surffile,nsurf,surf)
      call pdbread(pdbfile,mode,natom,atomch,resch,coord)
c
      write(*,'(A,A50)')'PDB file : ', pdbfile
      write(*,'(A,I5)') 'Atoms    : ', natom
      write(*,'(A,I5)') 'Surf pts : ', nsurf 
      if (eset.eq.1) write(*,'(A,F5.1)') 'TEMP     : ', temp
      write(*,'(A)') ' '
      write(*,'(A3,6A12)') 'SET', 'VdW area','VdW vol','SAS area', 
     &     'SAS vol ','SOL energy','Time '
      write(*,'(2A)') '---------------------------------------',
     &     '---------------------------------------'
      do set = 1, eset
c     
c     set = 1 - OONS
c     set = 2 - JRF_
c     set = 3 - WE92
c     set = 4 - SCH3
c     set = 5 - SCH4
c     
         call paraset(natom,atomch,resch,radi,
     &        gpara,hpara,cpara,vradi,rmax,debug,set) 
c     
         call sort(natom,coord,index,rank)
         call nlr(natom,coord,index,rmax-1.4,nls,nll)
         call nsol(natom,nls,nll,index,rank,coord,vradi,nsurf,
     &        surf,gpara,hpara,cpara,area,volume,esol,temp)
         varea = area
         vvolume = volume
         call nlr(natom,coord,index,rmax,nls,nll)
         time1 = etime(tar1)
         call nsol(natom,nls,nll,index,rank,coord,radi,nsurf,
     &        surf,gpara,hpara,cpara,area,volume,esol,temp)
         time2 = etime(tar2)
         write(*,'(I2,6F12.3)') set, varea, vvolume, area, volume,
     &           esol, time2 - time1
      end do
      write(*,'(2A)') '---------------------------------------',
     &     '---------------------------------------'
      write(*,'(A)') ' '
      write(*,'(A)') 'SET 1: [OONS] Ooi et al., 1987'
      write(*,'(A)') '    2: [JRF_] Vila et al., 1991'
      write(*,'(A)') '    3: [WE92] Wesson and Eisenberg, 1992'
      write(*,'(A)') '    4: [SCH3] Wesson and Eisenberg, 1992'
      write(*,'(A)') '    5: [SCH4] Schiffer et al., 1993'
c     
      stop
      end
c     
      subroutine nsol(natom,nls,nll,index,rank,coord,
     &     radi,nsurf,surf,gpara,hpara,cpara,area,volume,esol,temp)
c     
      parameter(maxatom=10000)
      parameter(maxsurf=10000)
      integer natom, atom, atom1, nsurf, count, l, i
      integer nls(maxatom,3), nll(maxatom,3), index(maxatom,3)
      integer rank(maxatom,3),key
      integer six
      double precision coord(maxatom,3), surf(maxsurf,3)
      double precision radi(maxatom), atarea(maxatom), radi2(maxatom),
     &     atvol(maxatom)
      double precision dx, dy, dz, dd, rr, sdx, sdy, sdz, sdd
      double precision gpara(maxatom), hpara(maxatom), cpara(maxatom)
      double precision gh, hh, ch, area, volume, esol, temp, temp0
      logical surfc(maxsurf)
      double precision xave, xmin, xmax, yave, ymin, ymax, zave, zmin,
     &     zmax
c     
      xave = coord(1,1)
      yave = coord(1,2)
      zave = coord(1,3)
      do i = 2, natom
         xave = xave + coord(i,1)
         yave = yave + coord(i,2)
         zave = zave + coord(i,3)
      end do
      xave = xave / dble(natom)
      yave = yave / dble(natom)
      zave = zave / dble(natom)

      do i = 1, natom
         radi2(i) = radi(i) * radi(i)
      end do
c
      area = 0.0d0
      volume = 0.0d0
      esol = 0.0d0
      gh = 0.0d0
      hh = 0.0d0
      ch = 0.0d0
      if(temp.eq.0) temp = 0.0000000000001d0
      temp0 = 298.0d0
      key = 1
      do atom = 1, natom
         six = rank(atom,key)
         do l = 1, nsurf
            surfc(l) = .true.
         end do
         do i = nls(six,key), nll(six,key)
            atom1 = index(i,key)
            if(atom.ne.atom1) then
               dx = coord(atom,1) - coord(atom1,1)
               dy = coord(atom,2) - coord(atom1,2)
               dz = coord(atom,3) - coord(atom1,3)
               rr = radi(atom) + radi(atom1)
               dd = dx * dx + dy * dy + dz * dz
               rr = rr * rr
               if (dd.le.rr) then
                  count = 0
                  do l = 1, nsurf
                     if(surfc(l)) then
                        count = count + 1
                        sdx = dx + surf(l,1) * radi(atom) 
                        sdy = dy + surf(l,2) * radi(atom) 
                        sdz = dz + surf(l,3) * radi(atom) 
                        sdd = sdx * sdx + sdy * sdy + sdz * sdz
                        rr = radi2(atom1)
                        if (sdd.le.rr) then
                           surfc(l) = .false.
                        end if
                     end if
                  end do
                  if(count.eq.0) goto 120
               end if
            end if
         end do
         count = 0
         dx = 0.0d0
         dy = 0.0d0
         dz = 0.0d0
         do l = 1, nsurf
            if(surfc(l)) then
               count = count + 1
               dx = dx + surf(l,1)
               dy = dy + surf(l,2)
               dz = dz + surf(l,3)
            end if
         end do
 120     continue
         sdr =  4.0d0 * 3.1415926535d0 
     &        * radi(atom) * radi(atom)
     &        / dble(nsurf)
         atarea(atom) = sdr * dble(count)
         atvol(atom) = sdr / 3.0d0
     &        * ((coord(atom,1) - xave)*dx
     &        + (coord(atom,2) - yave)*dy
     &        + (coord(atom,3) - zave)*dz
     &        +  dble(count) * radi(atom))
         if(count.eq.0) atvol(atom) = 0.0d0
         area = area + atarea(atom)
         volume = volume + atvol(atom)
         gh = gh + atarea(atom) * gpara(atom)
         hh = hh + atarea(atom) * hpara(atom)
         ch = ch + atarea(atom) * cpara(atom)
      end do
      esol =  (temp/temp0) * gh + hh * (1 - temp/temp0)
     &        - ch * (temp * dlog(temp/temp0) + temp0 - temp)
      return
      end
c     
      subroutine nlr(natom,coord,index,rmax,nls,nll)
c     
c     Neighbor list reduction by x-,y-,z- distance of each atoms.
c     Using sorted index, reduction can be shortly finished.
c     
      parameter(maxatom=10000)
      integer natom, si, sj, i, k
      integer nls(maxatom,3), nll(maxatom,3), index(maxatom,3)
      double precision coord(maxatom,3), x0, x1, d, rmax
      double precision rmax2
c     
      rmax2 = rmax * 2.0d0
c     
c     For each atom index(si,k) to the x-,y-,z-coordinates,
c     the farest atoms nls(si,k) and nll(si,k) are extracted.
c     
      do k = 1, 3
         do i = 1, natom
            nls(i,k) = 0
            nll(i,k) = 0
         end do
c     
         do si = 1, natom
            x0 = coord(index(si,k),k)
            do sj = si + 1, natom
               x1 = coord(index(sj,k),k)
c     x1 is always larger than x0.
               d = x1 - x0
               if(d.le.rmax2) then
                  nll(si,k) = sj
               else
                  goto 20
               end if
            end do
 20         continue
            do sj = si - 1, 1, -1
               x1 = coord(index(sj,k),k)
               d = x0 - x1
c     x1 is always smaller than x0.
               if(d.le.rmax2) then
                  nls(si,k) = sj
               else
                  goto 10
               end if
            end do
 10         continue
         end do
         nls(1,k) = 2
         nll(natom,k) = natom - 1
      end do
c     
      return
      end
c     
      subroutine sort(natom,coord,index,rank)
c     
c     index(#,[x,y,z]) is sorted atom no. by x-,y-,z-coordinates
c     for example, index(1,1) indicates an atom that has x-min,
c     index(natom,2) indicates an atom that has y-max.
c
c     rank(n,x) is rank of atom n by x
c
      parameter(maxatom=10000)
      integer i, j, k, tempi
      integer natom, index(maxatom,3), rank(maxatom,3)
      double precision coord(maxatom,3), a
c     
      do k = 1, 3
         do i = 1, natom
            index(i,k) = i
         end do
         do i = 2, natom
            tempi = index(i,k)
            a = coord(tempi,k)
            do j = i - 1, 1, -1
               if(coord(index(j,k),k).le.a) goto 10
               index(j+1,k) = index(j,k)
            end do
            j = 0
 10         index(j+1,k) = tempi
         end do
         do i = 1, natom
            rank(index(i,k),k) = i
         end do
      end do
c     
      return
      end
c     
      subroutine pdbread(pdbfile,mode,natom,atomch,resch,coord)
c     
c     pdb file read subroutine
c     
c     mode: 1 - all atoms
c     mode: 2 - all atoms except hydrogen atoms
c     mode: 3 - main chain only
c     
      parameter(maxatom=10000)
      integer mode, atomno(maxatom), resno(maxatom), natom, i
      character subunit(maxatom)*1, resch(maxatom)*3, atomch(maxatom)*4,
     &     u*1 ,unk1*1 ,unk2*1, tre*3, pdbfile*512, line*80
      double precision coord(maxatom,3)
c     
      open(unit=81,file=pdbfile,form='formatted',status='old',err=100)
      i = 0
 200  read(81,'(a)',end=300)line
      if(line(1:6).eq.'ATOM ') then
         i = i + 1
         read (line(7:80),'(i5,1x,4a,1x,a,i4,a,3x,3f8.3,1x)',end=300) 
     &        atomno(i),u,tre,unk1,resch(i),subunit(i),resno(i),
     &        unk2,coord(i,1),coord(i,2),coord(i,3)
         if ((u.eq.' ').or.(u.eq.'1').or.(u.eq.'2').or.(u.eq.'3')) then
            atomch(i) = tre//u
         else
            atomch(i) = u//tre
         end if
         if ((unk1.ne.' ').or.(unk2.ne.' ')) then
            write(*,'(3a)')'!pdb data error',unk1,unk2
            write(*,'(a)') line
         end if
         if(mode.eq.2) then
            if(atomch(i)(1:1).eq.'H') then
               i = i - 1
            end if
         end if
         if(mode.eq.3) then
            if((atomch(i)(1:2).ne.'N ')
     &           .and.(atomch(i)(1:2).ne.'CA')
     &           .and.(atomch(i)(1:2).ne.'C ')
     &           .and.(atomch(i)(1:2).ne.'O ')) then
               i = i - 1
            end if
         end if
      end if
      go to 200
 300  continue
      close(81)
      natom=i
      return
 100  write(*,'(a,a)') 'Cannot open ', pdbfile
      stop
      end
c     
      subroutine surfread(surffile,nsurf,surf)
c
      parameter(maxsurf=10000)
      integer nsurf, i
      double precision surf(maxsurf,3)
      character surffile*20, line*80
c
      open(unit=82,file=surffile,form='formatted',status='old',err=100)
      go to 200
 100  write(*,'(a,a)') surffile,' cannot open.'
      stop
 200  read(82,'(a)')line
      if (line(1:1).eq.'#')then
         go to 200
      end if
      i = 1
      read(line(1:80),'(3f20.10)') surf(i,1), surf(i,2), 
     &     surf(i,3)
 300  i = i + 1
      read(82,'(3f20.10)',end=400) surf(i,1), surf(i,2), 
     &     surf(i,3)
      go to 300
 400  continue
      nsurf = i - 1
      close(82)
      return
      end
c
      subroutine paraset(natom,atomch,resch,radi,gpara,hpara,cpara,
     &                   vradi,rmax,debug,set)
c     
      parameter(maxatom=10000)
      integer natom, i, set
      logical debug
      character resch(maxatom)*3, atomch(maxatom)*4, setname*4
      double precision radi(maxatom), gpara(maxatom), vradi(maxatom)
      double precision hpara(maxatom), cpara(maxatom), rmax
C
C                          radi    OONS     JRF_   OONS(H)  OONS(C/1000)
C     Aliphatic         C  2.00    0.008    0.216   -0.026  0.370
C     Aromatic          C  1.75   -0.008   -0.678   -0.038  0.296
C     Carbonyl/carboxyl C  1.55    0.427   -0.732    0.413  0.613
C     Amide             N  1.55   -0.132   -0.312   -0.192 -0.012
C     Carbonyl/carboxyl O  1.40   -0.038   -0.262   -0.032 -0.228
C     Hydroxyl          O  1.40   -0.172   -0.910   -0.238  0.008
C     Sulfer            S  2.00   -0.021   -0.281   -0.031 -0.001
C
C                          radi    WE92     SCH3     SCH4
C     Carbon            C  1.90    0.012    0.004    0.0325
C     Non-charged       O  1.40   -0.116   -0.113   -0.0095
C     Charged           O  1.40   -0.175   -0.116   -0.280
C     Non-charged       N  1.70   -0.116   -0.113   -0.0175
C     Charged           N  1.70   -0.186   -0.169   -0.2175
C     Sulfer            S  1.80   -0.018   -0.017   -0.009
C
C     OONS: Ooi et al., 1987
C     JRF_: Vila et al., 1991
C     WE92: Wesson and Eisenberg, 1992
C     SCH3: Wesson and Eisenberg, 1992
C     SCH4: Schiffer et al., 1993
C
      double precision pc1(5), pc2(5), pc3(5), po1(5), po2(5), pn1(5), 
     &     pn2(5), ps0(5), rc1(5), rc2(5), rc3(5), ro1(5), ro2(5),
     &     rn1(5), rn2(5), rs0(5), rwa, ph0, rh0,
     &     qc1(2), qc2(2), qc3(2), qn1(2), qn2(2), qo1(2), qo2(2),
     &     qs0(2)
C
C     SET = 1 : OONS
C     SET = 2 : JRF_
C     SET = 3 : WE92
C     SET = 4 : SCH3
C     SET = 5 : SCH4
C
      data ph0/0.0d0/
      data pc1/ 0.008d0, 0.216d0, 0.012d0, 0.004d0, 0.0325d0/
      data pc2/-0.008d0,-0.678d0, 0.012d0, 0.004d0, 0.0325d0/
      data pc3/ 0.427d0,-0.732d0, 0.012d0, 0.004d0, 0.0325d0/
      data pn1/-0.132d0,-0.312d0,-0.116d0,-0.113d0,-0.0175d0/
      data pn2/-0.132d0,-0.312d0,-0.186d0,-0.169d0,-0.2175d0/
      data po1/-0.038d0,-0.262d0,-0.116d0,-0.113d0,-0.0095d0/
      data po2/-0.172d0,-0.910d0,-0.175d0,-0.116d0,-0.2800d0/
      data ps0/-0.021d0,-0.281d0,-0.018d0,-0.017d0,-0.0090d0/
c
      data qc1/-0.026d0, 0.000370d0/
      data qc2/-0.038d0, 0.000296d0/
      data qc3/ 0.413d0, 0.000613d0/
      data qn1/-0.192d0,-0.000012d0/
      data qn2/-0.192d0,-0.000012d0/
      data qo1/-0.032d0,-0.000228d0/
      data qo2/-0.238d0, 0.000008d0/
      data qs0/-0.031d0,-0.000001d0/
c
      data rwa/1.40d0/
      data rh0/0.00d0/
      data rc1/2.00d0,2.00d0,1.90d0,1.90d0,1.90d0/
      data rc2/1.75d0,1.75d0,1.90d0,1.90d0,1.90d0/
      data rc3/1.55d0,1.55d0,1.90d0,1.90d0,1.90d0/
      data ro1/1.40d0,1.40d0,1.40d0,1.40d0,1.40d0/
      data ro2/1.40d0,1.40d0,1.40d0,1.40d0,1.40d0/
      data rn1/1.55d0,1.55d0,1.70d0,1.70d0,1.70d0/
      data rn2/1.55d0,1.55d0,1.70d0,1.70d0,1.70d0/
      data rs0/2.00d0,2.00d0,1.80d0,1.80d0,1.80d0/
c
      do i=1, natom
        radi(i) = 0.0d0
        gpara(i) = 0.0d0
        hpara(i) = 0.0d0
        cpara(i) = 0.0d0
        vradi(i) = 0.0d0
      end do
c
      if(set.le.2) then
          rmax = 3.4d0
      else
          rmax = 3.3d0
      end if
c
      if (debug) then
         open(unit=30,file='nsol.out',form='formatted',status='unknown')
         if (set.eq.1) then
            setname = 'oons'
         else if (set.eq.2) then
            setname = 'jrf_'
         else if (set.eq.3) then 
            setname = 'we92'
         else if (set.eq.4) then
            setname = 'sch3'
         else if (set.eq.5) then
            setname = 'sch4'
         else
            setname = 'unk '
         end if
         write(30,'(a,a)') 'parameter - ',setname
         write(30,'(a3,a,a4,a3,3a8)')
     &        'no ',' ','atm ','res','radi ','vradi','para '
      end if
c
      do i = 1, natom
c     hydrogen atom
         if(atomch(i)(1:1).eq.'H') then
            gpara(i) = ph0
            hpara(i) = ph0
            cpara(i) = ph0
            radi(i) = rh0
c     sulfer
         else if (atomch(i)(1:1).eq.'S') then
            gpara(i) = ps0(set)
            hpara(i) = qs0(1)
            cpara(i) = qs0(2)
            radi(i) = rs0(set)
c     nitrogen
         else if (atomch(i)(1:1).eq.'N') then
c     charged nitrogen
            if (
     &           ((atomch(i)(2:3).eq.'H2')
     &           .and.(resch(i)(1:3).eq.'ARG'))
     &           .or.
     &           (atomch(i)(2:2).eq.'Z')
     &           .or.
     &           ((atomch(i)(2:3).eq.'E2')
     &           .and.(resch(i)(1:2).eq.'HI'))
     &           ) then
               gpara(i) = pn2(set)
               hpara(i) = qn2(1)
               cpara(i) = qn2(2)
               radi(i) = rn2(set)
c     non-charged nitrogen
            else
               gpara(i) = pn1(set)
               hpara(i) = qn1(1)
               cpara(i) = qn1(2)
               radi(i) = rn1(set)
            end if
c     oxgen
         else if (atomch(i)(1:1).eq.'O') then
            if(set.ge.3) then
c     we92, sch3, sch4
c     charged oxgen
               if(
     &              ((resch(i)(1:3).eq.'ASP')
     &              .and.(atomch(i)(1:3).eq.'OD2'))
     &              .or.
     &              ((resch(i)(1:3).eq.'GLU')
     &              .and.(atomch(i)(1:3).eq.'OE2'))
     &              ) then
                  gpara(i) = po2(set)
                  hpara(i) = qo2(1)
                  cpara(i) = qo2(2)
                  radi(i) = ro2(set)
c     non-charged oxgen
               else if (atomch(i)(1:3).eq.'OXT') then
                  gpara(i) = po1(set)
                  hpara(i) = qo1(1)
                  cpara(i) = qo1(2)
                  radi(i) = ro1(set)
                else
                  gpara(i) = po1(set)
                  hpara(i) = qo1(1)
                  cpara(i) = qo1(2)
                  radi(i) = ro1(set)
               end if
            else
c     oons, jrf
c     carbonyl/carboxyl oxgen
               if (atomch(i)(2:2).eq.' ') then
                  gpara(i) = po1(set)
                  hpara(i) = qo1(1)
                  cpara(i) = qo1(2)
                  radi(i) = ro1(set)
               else if (atomch(i)(1:3).eq.'OXT') then
                  gpara(i) = po1(set)
                  hpara(i) = qo1(1)
                  cpara(i) = qo1(2)
                  radi(i) = ro1(set)
               else if((atomch(i)(2:3).eq.'E1')
     &              .or.
     &              (atomch(i)(2:3).eq.'D1')
     &              ) then
                  gpara(i) = po1(set)
                  hpara(i) = qo1(1)
                  cpara(i) = qo1(2)
                  radi(i) = ro1(set)
               else
                  gpara(i) = po2(set)
                  hpara(i) = qo2(1)
                  cpara(i) = qo2(2)
                  radi(i) = ro2(set)
               end if
            end if
c     oxgen end
c     carbon
         else if (atomch(i)(1:1).eq.'C') then
            if (atomch(i)(2:2).eq.' ') then
               gpara(i) = pc3(set)
               hpara(i) = qc3(1)
               cpara(i) = qc3(2)
               radi(i) = rc3(set)
            else if ((atomch(i)(2:2).eq.'A')
     &              .or.(atomch(i)(2:2).eq.'B')) then
               gpara(i) = pc1(set)
               hpara(i) = qc1(1)
               cpara(i) = qc1(2)
               radi(i) = rc1(set)
            else if ((atomch(i)(2:2).eq.'H')
     &              .or.(atomch(i)(2:2).eq.'Z')) then
               gpara(i) = pc2(set)
               hpara(i) = qc2(1)
               cpara(i) = qc2(2)
               radi(i) = rc2(set)
            else if (atomch(i)(2:2).eq.'D') then
               if (resch(i)(1:2).eq.'GL') then
                  gpara(i) = pc3(set)
                  hpara(i) = qc3(1)
                  cpara(i) = qc3(2)
                  radi(i) = rc3(set)
               else if ((resch(i)(1:3).eq.'ARG')
     &              .or.(resch(i)(1:3).eq.'ILE')
     &              .or.(resch(i)(1:3).eq.'LYS')
     &              .or.(resch(i)(1:3).eq.'PRO')
     &              .or.(resch(i)(1:3).eq.'LEU')) then
                  gpara(i) = pc1(set)
                  hpara(i) = qc1(1)
                  cpara(i) = qc1(2)
                  radi(i) = rc1(set)
               else 
                  gpara(i) = pc2(set)
                  hpara(i) = qc2(1)
                  cpara(i) = qc2(2)
                  radi(i) = rc2(set)
               end if
            else if (atomch(i)(2:2).eq.'E') then
               if (atomch(i)(3:3).eq.' ') then
                  gpara(i) = pc1(set)
                  hpara(i) = qc1(1)
                  cpara(i) = qc1(2)
                  radi(i) = rc1(set)
               else 
                  gpara(i) = pc2(set)
                  hpara(i) = qc2(1)
                  cpara(i) = qc2(2)
                  radi(i) = rc2(set)
               end if
            else if (atomch(i)(2:2).eq.'G') then
               if (resch(i)(1:2).eq.'AS') then
                  gpara(i) = pc3(set)
                  hpara(i) = qc3(1)
                  cpara(i) = qc3(2)
                  radi(i) = rc3(set)
               else if ((resch(i)(1:2).eq.'HI')
     &                 .or.(resch(i)(1:3).eq.'PHE')
     &                 .or.(resch(i)(1:3).eq.'TRP')
     &                 .or.(resch(i)(1:3).eq.'TYR')) then
                  gpara(i) = pc2(set)
                  hpara(i) = qc2(1)
                  cpara(i) = qc2(2)
                  radi(i) = rc2(set)
               else 
                  gpara(i) = pc1(set)
                  hpara(i) = qc1(1)
                  cpara(i) = qc1(2)
                  radi(i) = rc1(set)
               end if
            end if
         else 
            gpara(i) = 0.0d0
            hpara(i) = 0.0d0
            cpara(i) = 0.0d0
            radi(i) = 0.0d0
            write(*,'(i3,1x,a4,a3,a)') i, atomch(i), resch(i),
     &           ' is not in parameter database. set to 0.0d0.'
         end if
         vradi(i) = radi(i)
         if (radi(i).ne.0.0d0) radi(i) = radi(i) + rwa
         if (debug) then
            write(30,'(i3,a,a4,a3,2f8.3,3f8.4)') 
     &           i,' ',atomch(i),resch(i),radi(i),vradi(i),gpara(i),
     &           hpara(i), cpara(i)
         end if
      end do
      if (debug) close(30)
      return
      end