Changeset c076b43 for enylun.f


Ignore:
Timestamp:
09/24/08 15:27:03 (16 years ago)
Author:
sandipan <sandipan@…>
Branches:
master
Children:
4fd4338
Parents:
a52fb83
Message:

Updated version of Lund force field functions

The Lund force field functions were debugged. A new library file
lib.lun was added to describe the geometrical parameters of
protein molecules consistent with the force field. With this
library file defining the bond lengths and bond angles, energy
values calculated with SMMP with the Lund force field can be
directly compared with those calculated with PROFASI.

git-svn-id: ​svn+ssh://svn.berlios.de/svnroot/repos/smmp/trunk@17 26dc1dd8-5c4e-0410-9ffe-d298b4865968

File:
1 edited

Legend:

Unmodified
Added
Removed

  • enylun.f

    ra52fb83 rc076b43  
    1010!                      Jan H. Meinke, Sandipan Mohanty
    1111!
     12      subroutine set_local_constr(ires,frst,lst,root)
     13      include 'INCL.H'
     14      include 'incl_lund.h'
     15          integer ires,frst,lst,root
     16          do iat1=iCa(ires)+frst,iCa(ires)+lst
     17             do iat2=iat1+1,iCa(ires)+lst
     18                matcon(iat2-iat1,iat1)=0
     19                matcon(iat1-iat2,iat2)=0
     20             enddo
     21          enddo
     22          iat1=iCa(ires)+root
     23          do iat2=iCa(ires)+root,iCa(ires)+lst
     24                matcon(iat2-iat1,iat1)=0
     25                matcon(iat1-iat2,iat2)=0
     26          enddo
     27      end subroutine set_local_constr
     28
    1229      subroutine init_lundff
    1330      include 'INCL.H'
    … …  
    1936!     Some parameters in the Lund force field.
    2037!     The correspondence between internal energy scale and kcal/mol
    21       eunit=1.3315
     38      eunit=1.3315d0
     39!      eunit=1.0
    2240!     Bias
    23       kbias=100.0*eunit
     41      kbias=100.0d0*eunit
    2442!      print *,'Bias'
    2543!     Hydrogen bonds
    26       epshb1=3.1*eunit
    27       epshb2=2.0*eunit
    28       sighb=2.0
    29       cthb=4.5
     44      epshb1=3.1d0*eunit
     45      epshb2=2.0d0*eunit
     46      sighb=2.0d0
     47      cthb=4.5d0
    3048      cthb2=cthb*cthb
    31       powa=0.5
    32       powb=0.5
    33       blhb=-30.0*(((sighb/cthb)**10-(sighb/cthb)**12))/cthb2
     49      powa=0.5d0
     50      powb=0.5d0
     51      blhb=-30.0d0*(((sighb/cthb)**10-(sighb/cthb)**12))/cthb2
    3452      alhb=-(5*((sighb/cthb)**12)-6*((sighb/cthb)**10))-blhb*cthb2
    3553      sighb2=sighb*sighb
    36       cdon=1.0
    37       cacc=(1.0/1.23)**powb
    38       csacc=(1.0/1.25)**powb
     54      cdon=1.0d0
     55      cacc=(1.0d0/1.23d0)**powb
     56      csacc=(1.0d0/1.25d0)**powb
    3957!      print *,'Hydrogen bonds'
    4058!     Hydrophobicity
    4159!      print *,'Hydrophobicity with nhptyp = ',nhptyp
    4260
    43       hpstrg(1)=0.0*eunit
    44       hpstrg(2)=0.1*eunit
    45       hpstrg(3)=0.1*eunit
    46       hpstrg(4)=0.1*eunit
    47       hpstrg(5)=0.9*eunit
    48       hpstrg(6)=2.8*eunit
    49       hpstrg(7)=0.1*eunit
    50       hpstrg(8)=2.8*eunit
    51       hpstrg(9)=3.2*eunit
     61      hpstrg(1)=0.0d0*eunit
     62      hpstrg(2)=0.1d0*eunit
     63      hpstrg(3)=0.1d0*eunit
     64      hpstrg(4)=0.1d0*eunit
     65      hpstrg(5)=0.9d0*eunit
     66      hpstrg(6)=2.8d0*eunit
     67      hpstrg(7)=0.1d0*eunit
     68      hpstrg(8)=2.8d0*eunit
     69      hpstrg(9)=3.2d0*eunit
    5270
    5371      do i=1,mxrs
    … …  
    6482     &           .or.(mynm.eq.'pron').or.(mynm.eq.'pro+')) then
    6583            prlvr=.true.        ! residue i is a proline variant
     84            print *, 'proline variant ',mynm,i
    6685         else
    6786            prlvr=.false.
    … …  
    86105            ihpat(i,2)=iCa(i)+5
    87106            ihpat(i,3)=iCa(i)+7
    88             ihpat(i,3)=iCa(i)+11
     107            ihpat(i,4)=iCa(i)+11
    89108         else if (mynm.eq.'ile') then
    90109            nhpat(i)=4
    … …  
    92111            ihpat(i,2)=iCa(i)+4
    93112            ihpat(i,3)=iCa(i)+8
    94             ihpat(i,3)=iCa(i)+11
     113            ihpat(i,4)=iCa(i)+11
    95114         else if (mynm.eq.'met') then
    96115            nhpat(i)=4
    … …  
    98117            ihpat(i,2)=iCa(i)+5
    99118            ihpat(i,3)=iCa(i)+8
    100             ihpat(i,3)=iCa(i)+9
     119            ihpat(i,4)=iCa(i)+9
    101120         else if (mynm.eq.'phe') then
    102121            nhpat(i)=6
    … …  
    104123            ihpat(i,2)=iCa(i)+6
    105124            ihpat(i,3)=iCa(i)+8
    106             ihpat(i,3)=iCa(i)+10
    107             ihpat(i,3)=iCa(i)+12
    108             ihpat(i,3)=iCa(i)+14
     125            ihpat(i,4)=iCa(i)+10
     126            ihpat(i,5)=iCa(i)+12
     127            ihpat(i,6)=iCa(i)+14
    109128         else if (mynm.eq.'tyr') then
    110129            nhpat(i)=6
    … …  
    112131            ihpat(i,2)=iCa(i)+6
    113132            ihpat(i,3)=iCa(i)+8
    114             ihpat(i,3)=iCa(i)+10
    115             ihpat(i,3)=iCa(i)+13
    116             ihpat(i,3)=iCa(i)+15
     133            ihpat(i,4)=iCa(i)+10
     134            ihpat(i,5)=iCa(i)+13
     135            ihpat(i,6)=iCa(i)+15
    117136         else if (mynm.eq.'trp') then
    118137            nhpat(i)=6
    … …  
    120139            ihpat(i,2)=iCa(i)+11
    121140            ihpat(i,3)=iCa(i)+13
    122             ihpat(i,3)=iCa(i)+15
    123             ihpat(i,3)=iCa(i)+17
    124             ihpat(i,3)=iCa(i)+19
     141            ihpat(i,4)=iCa(i)+15
     142            ihpat(i,5)=iCa(i)+17
     143            ihpat(i,6)=iCa(i)+19
    125144         endif
    126145      enddo
    … …  
    132151      exvcut2=exvcut*exvcut
    133152
    134       sigsa(1)=1.0  ! hydrogen
    135       sigsa(2)=1.0  ! hydrogen
    136       sigsa(3)=1.0  ! hydrogen
    137       sigsa(4)=1.0  ! hydrogen
    138       sigsa(5)=1.0  ! hydrogen
    139       sigsa(6)=1.0  ! hydrogen
    140       sigsa(7)=1.75 ! carbon
    141       sigsa(8)=1.75 ! carbon
    142       sigsa(9)=1.75 ! carbon
    143       sigsa(10)=1.42 ! oxygen
    144       sigsa(11)=1.42 ! oxygen
    145       sigsa(12)=1.42 ! oxygen
    146       sigsa(13)=1.55 ! nitrogen
    147       sigsa(14)=1.55 ! nitrogen
    148       sigsa(15)=1.55 ! nitrogen
    149       sigsa(16)=1.77 ! sulfur
    150       sigsa(17)=1.0  ! hydrogen
    151       sigsa(18)=1.75 ! carbon
     153      sigsa(1)=1.0d0  ! hydrogen
     154      sigsa(2)=1.0d0  ! hydrogen
     155      sigsa(3)=1.0d0  ! hydrogen
     156      sigsa(4)=1.0d0  ! hydrogen
     157      sigsa(5)=1.0d0  ! hydrogen
     158      sigsa(6)=1.0d0  ! hydrogen
     159      sigsa(7)=1.75d0 ! carbon
     160      sigsa(8)=1.75d0 ! carbon
     161      sigsa(9)=1.75d0 ! carbon
     162      sigsa(10)=1.42d0 ! oxygen
     163      sigsa(11)=1.42d0 ! oxygen
     164      sigsa(12)=1.42d0 ! oxygen
     165      sigsa(13)=1.55d0 ! nitrogen
     166      sigsa(14)=1.55d0 ! nitrogen
     167      sigsa(15)=1.55d0 ! nitrogen
     168      sigsa(16)=1.77d0 ! sulfur
     169      sigsa(17)=1.0d0  ! hydrogen
     170      sigsa(18)=1.75d0 ! carbon
    152171
    153172      do i=1,mxtyat
    154173         do j=1,mxtyat
    155174            sig2lcp(i,j)=(sigsa(i)+sigsa(j))**2
    156             asalcp(i,j)=-7*((sig2lcp(i,j)/exvcut2)**6.0)
    157             bsalcp(i,j)=6*((sig2lcp(i,j)/exvcut2)**6.0)/exvcut2
     175            asalcp(i,j)=-7*((sig2lcp(i,j)/exvcut2)**6.0d0)
     176            bsalcp(i,j)=6*((sig2lcp(i,j)/exvcut2)**6.0d0)/exvcut2
    158177         enddo
    159178      enddo
    160179!      print *,'Local pair excluded volume constants'
    161180
    162       exvlam=0.75
     181      exvlam=0.75d0
    163182      exvcutg=exvcut*exvlam
    164183      exvcutg2=exvcutg*exvcutg
    … …  
    233252            matcon(iat1-iCa(irsml1(iml)),iCa(irsml1(iml)))=0
    234253            matcon(iCa(irsml1(iml))-iat1,iat1)=0
    235 
     254            do iat2=iat1,iCa(irsml1(iml))-1
     255                matcon(iat2-iat1,iat1)=0
     256                matcon(iat1-iat2,iat2)=0
     257            enddo
    236258            matcon(iat1-iCa(irsml1(iml))-1,iCa(irsml1(iml))+1)=2
    237259            matcon(iCa(irsml1(iml))+1-iat1,iat1)=2
    … …  
    259281            mynm=seq(irs)
    260282            call tolost(mynm)
    261             if ((mynm.eq.'pro').or.(mynm.eq.'cpro')
    262      &              .or.(mynm.eq.'cpru').or.(mynm.eq.'prou')
    263      &              .or.(mynm.eq.'pron').or.(mynm.eq.'pro+')) then
     283            if ((mynm.eq.'pro')) then
    264284               prlvr=.true.     ! residue i is a proline variant
    265285            else
    266286               prlvr=.false.
    267287            endif
    268 
    269             if ((mynm.eq.'arg').or.(mynm.eq.'arg+')) then
    270                do iat1=iatoff+iatrs1(irs)+13,iatrs2(irs)-2-iatmrg
    271                   do iat2=iat1+1,iatrs2(irs)-2-iatmrg
    272                      matcon(iat2-iat1,iat1)=0
    273                      matcon(iat1-iat2,iat2)=0
    274                   enddo
    275                enddo
    276             else if ((mynm.eq.'his').or.(mynm.eq.'hise')
    277      &              .or.(mynm.eq.'hisd').or.(mynm.eq.'his+')) then
    278                do iat1=iatoff+iatrs1(irs)+7,iatrs2(irs)-2-iatmrg
    279                   do iat2=iat1+1,iatrs2(irs)-2-iatmrg
    280                      matcon(iat2-iat1,iat1)=0
    281                      matcon(iat1-iat2,iat2)=0
    282                   enddo
    283                enddo
     288            if ((mynm.eq.'asn')) then
     289                call set_local_constr(irs,5,9,2)
     290            else if ((mynm.eq.'gln')) then
     291                call set_local_constr(irs,8,12,5)
     292            else if ((mynm.eq.'arg')) then
     293                call set_local_constr(irs,11,19,8)
     294            else if ((mynm.eq.'his')) then
     295                call set_local_constr(irs,5,12,2)
    284296            else if (mynm.eq.'phe') then
    285                do iat1=iatoff+iatrs1(irs)+7,iatrs2(irs)-2-iatmrg
    286                   do iat2=iat1+1,iatrs2(irs)-2-iatmrg
    287                      matcon(iat2-iat1,iat1)=0
    288                      matcon(iat1-iat2,iat2)=0
    289                   enddo
    290                enddo
     297                call set_local_constr(irs,5,15,2)
    291298            else if (mynm.eq.'tyr') then
    292                do iat1=iatoff+iatrs1(irs)+7,iatrs2(irs)-2-iatmrg
    293                   do iat2=iat1+1,iatrs2(irs)-2-iatmrg
    294                      if (iat2.ne.iatrs1(irs)+15) then
    295                         matcon(iat2-iat1,iat1)=0
    296                         matcon(iat1-iat2,iat2)=0
    297                      endif
    298                   enddo
    299                enddo
     299                call set_local_constr(irs,5,16,2)
     300                iat1=iCa(irs)+12 ! H_h
     301                iat2=iCa(irs)+13 ! C_e2
     302                iat3=iCa(irs)+8  ! C_e1
     303                matcon(iat2-iat1,iat1)=2
     304                matcon(iat1-iat2,iat2)=2
     305                matcon(iat3-iat1,iat1)=2
     306                matcon(iat1-iat3,iat3)=2
    300307            else if (mynm.eq.'trp') then
    301                do iat1=iatoff+iatrs1(irs)+7,iatrs2(irs)-2-iatmrg
    302                   do iat2=iat1+1,iatrs2(irs)-2-iatmrg
    303                      matcon(iat2-iat1,iat1)=0
    304                      matcon(iat1-iat2,iat2)=0
    305                   enddo
    306                enddo
     308                call set_local_constr(irs,5,19,2)
    307309            else if (prlvr) then
    308310!           Proline. Many more distances are fixed because of the fixed
    309311!           phi angle
    310                do iat1=iatoff+iatrs1(irs),iatrs2(irs)-2-iatmrg
    311                   do iat2=iat1+1,iatrs2(irs)-2-iatmrg
    312                      matcon(iat2-iat1,iat1)=0
    313                      matcon(iat1-iat2,iat2)=0
    314                   enddo
    315                enddo
    316 !           distances to the C' atom of the previous residue are also fixed
    317                if (irs.ne.irsml1(iml)) then
    318                   iat1=iowat(iatrs1(irs))
    319                   do iat2=iatrs1(irs),iatrs2(irs)-1-iatmrg
    320                      matcon(iat2-iat1,iat1)=0
    321                      matcon(iat1-iat2,iat2)=0
    322                   enddo
    323                endif
     312                call set_local_constr(irs,-3,11,-3)
     313                matcon(iCa(irs-1)-iCa(irs)-8,iCa(irs)+8)=0
     314                matcon(iCa(irs)+8-iCa(irs-1),iCa(irs-1))=0
    324315            endif
    325316         enddo
     317
     318!        Distances fixed because of the constant omega angle
     319         do irs=irsml1(iml),irsml2(iml)-1
     320            matcon(iCa(irs+1)-iC(irs)-1,iC(irs)+1)=0 ! O_i -- Ca_{i+1}
     321            matcon(-iCa(irs+1)+iC(irs)+1,iCa(irs+1))=0 ! O_i -- Ca_{i+1}
     322
     323            matcon(iN(irs+1)-iC(irs),iC(irs)+1)=0    ! O_i -- H_{i+1}
     324            matcon(-iN(irs+1)+iC(irs),iN(irs+1)+1)=0    ! O_i -- H_{i+1}
     325
     326            matcon(iCa(irs+1)-iCa(irs),iCa(irs))=0   ! Ca_i -- Ca_{i+1}
     327            matcon(-iCa(irs+1)+iCa(irs),iCa(irs+1))=0   ! Ca_i -- Ca_{i+1}
     328
     329            matcon(iN(irs+1)+1-iCa(irs),iCa(irs))=0  ! Ca_i -- H_{i+1}
     330            matcon(-iN(irs+1)-1+iCa(irs),iN(irs+1)+1)=0  ! Ca_i -- H_{i+1}
     331         enddo
     332
    326333      enddo
    327334!     finished initializing matrix conmat
    328335!      print *,'Connections matrix'
    329 
    330336!     Local pair excluded volume
    331337      do i=1,mxml
    … …  
    354360            ilpst(iml+1)=ilp+1
    355361         endif
    356 !         print *,'molecule ',iml,' lc pair range ',ilpst(iml),ilpnd(iml)
    357 !         print *,'local pair list'
     362         print *,'molecule ',iml,' lc pair range ',ilpst(iml),ilpnd(iml)
     363!        print *,'local pair list'
    358364         do lci=ilpst(iml),ilpnd(iml)
    359365            iat1=lcp1(lci)
    360366            iat2=lcp2(lci)
    361 !            print *,lci,iat1,iat2,matcon(iat2-iat1,iat1)
     367!            print *,lci,iat1,iat2,ityat(iat1),ityat(iat2)
    362368         enddo
    363369      enddo
    … …  
    375381         ityp=1
    376382      else if ((mynm.eq.'val').or.(mynm.eq.'leu').or.(mynm.eq.'ile')
    377      &        .or.(mynm.eq.'met').or.(mynm.eq.'pro').or.(mynm.eq.'cpro')
    378      &        .or.(mynm.eq.'cpru').or.(mynm.eq.'prou')
    379      &        .or.(mynm.eq.'pron').or.(mynm.eq.'pro+')) then
     383     &        .or.(mynm.eq.'met').or.(mynm.eq.'pro')) then
    380384         ityp=2
    381385      else if ((mynm.eq.'phe').or.(mynm.eq.'tyr').or.(mynm.eq.'trp'))
    … …  
    444448      r4=r6*r6
    445449      r6=r6*r4
    446       evlu=((ca*ca/r2)**(0.5*powa))*((cb*cb/r2)**(0.5*powb))
     450      rdon2=(xat(d2)-xat(d1))**2+(yat(d2)-yat(d1))**2+
     451     & (zat(d2)-zat(d1))**2
     452      racc2=(xat(a2)-xat(a1))**2+(yat(a2)-yat(a1))**2+
     453     & (zat(a2)-zat(a1))**2
     454      evlu=((ca*ca/(r2*rdon2))**(0.5*powa))*
     455     & ((cb*cb/(r2*racc2))**(0.5*powb))
    447456      evlu=evlu*(r6*(5*r6-6*r4)+alhb+blhb*r2)
    448 !      print *,'found hbmm contribution ',evlu
     457!      print *,'found hbmm contribution ',i,j,epshb1,evlu
    449458      ehbmmrs=epshb1*evlu
    450459      return
    … …  
    495504         if (.not.prlvr) then 
    496505            do j=istres,indres
    497                if ((j.eq.istres).or.(j.eq.indres)) shbm2=0.5
    498                etmp=ehbmmrs(i,j)
    499                eyhb=eyhb+shbm1*shbm2*etmp
     506               if ((j.lt.(i-2)).or.(j.gt.(i+1))) then
     507                    shbm2=1.0
     508                    if ((j.eq.istres).or.(j.eq.indres)) shbm2=0.5
     509                    etmp=ehbmmrs(i,j)
     510                    eyhb=eyhb+shbm1*shbm2*etmp
     511                endif
    500512            enddo
    501513         endif
    … …  
    512524                  if (j.eq.(i+1)) etmp=0
    513525                  if (j.eq.(i+2)) etmp=0.5*etmp
     526!                  print *, 'hydrophobicity contribution ',etmp,i,j
    514527                  eysl=eysl+etmp
    515528               endif
    … …  
    541554            r6=r6*r6*r6
    542555            etmp1=(r6*r6+asalcp(iatt1,iatt2)+bsalcp(iatt1,iatt2)*r2)
    543             if (etmp1.ge.500) then
    544                print *,'local pair contribution ',iat1,iat2,etmp1
     556            if (etmp1.ge.0) then
     557!               print *,'local pair contribution ',iat1,iat2,iatt1,
     558!     & iatt2,sqrt(sig2lcp(iatt1,iatt2)),etmp1
    545559            endif
    546560            etmp=etmp+etmp1
    … …  
    592606         enddo
    593607      enddo
    594       sum=0
     608      ssum=0
    595609      do i=1,ni+nj
    596610         if (r2min(i).le.b2) then
    597611            if (r2min(i).lt.a2) then
    598                sum=sum+1
     612               ssum=ssum+1
    599613            else
    600                sum=sum+(b2-r2min(i))/(b2-a2)
     614               ssum=ssum+(b2-r2min(i))/(b2-a2)
    601615            endif
    602616         endif
    603       enddo
    604       ehp=-hpstrg((ihp1-1)*nhptyp+ihp2)*sum/(ni+nj)
     617!         hpstrgth=hpstrg((ihp1-1)*nhptyp+ihp2)
     618!         print *, 'hp diagnosis ',ni,nj,ssum/(ni+nj),hpstrgth
     619      enddo
     620      ehp=-hpstrg((ihp1-1)*nhptyp+ihp2)*ssum/(ni+nj)
    605621      return
    606622      end
    … …  
    617633      dimension isort(mxat),ngbr(mxat),locccl(mxat),incell(mxcell)
    618634      dimension icell(mxat)
     635      integer :: jx1, jy1, jz1
     636
    619637      if (nml.eq.0) then
    620638         istat=iatrs1(irsml1(1))
    … …  
    736754         lcell=locccl(icl)
    737755         ix=mod(lcell-1,ndx)
    738          iy=(mod(lcell-1,nxy)-ix)/ndx
    739          iz=(lcell-1-ix-ndx*iy)/nxy
    740 !         print *,'icl=',icl,'absolute index of cell = ',lcell
    741 !         print *,'iz,iy,ix = ',iz,iy,ix
    742 !     find all atoms in current cell and all its forward-going neighbours
    743          nex=min(ix+1,ndx-1)
    744          ney=min(iy+1,ndy-1)
    745          nez=min(iz+1,ndz-1)
     756         iz = (lcell - 1) / nxy
     757         iy = ((lcell - 1) - iz * nxy) / ndx
     758
     759c         print *,'icl=',icl,'absolute index of cell = ',lcell
     760c         print *,'iz,iy,ix = ',iz,iy,ix
     761c     find all atoms in current cell and all its forward-going neighbours
     762         nex=ix+1
     763         ney=iy+1
     764         nez=iz+1
    746765         nsame=0
    747766         nngbr=0
    748767         do jx=ix,nex
     768            if (jx.ge.ndx) then
     769               jx1=0
     770            else
     771               jx1=jx
     772            endif
    749773            do jy=iy,ney
     774               if (jy.ge.ndy) then
     775                  jy1=0
     776               else
     777                  jy1=jy
     778               endif
    750779               do jz=iz,nez
    751                   jcl=jx+ndx*jy+nxy*jz+1
     780                  if (jz.ge.ndz) then
     781                     jz1=0
     782                  else
     783                     jz1=jz
     784                  endif
     785                  jcl=jx1 + ndx*jy1 + nxy*jz1 + 1
     786!                  write(*,*)'jcl,jx1,jy1,jz1:', jcl,jx1,jy1,jz1
    752787                  do ii=incell(jcl)+1,incell(jcl+1)
    753 !                    count the total number of neighbours
     788!     count the total number of neighbours
    754789                     nngbr=nngbr+1
    755790                     if (jx.eq.ix.and.jy.eq.iy.and.jz.eq.iz) then
    756 !                    count how many neighbours are from the same cell
     791!     count how many neighbours are from the same cell
    757792                        nsame=nsame+1
    758793                     endif
    … …  
    766801!        1
    767802         jx=ix+1
     803         if (jx.ge.ndx) jx=0
    768804         jy=iy
    769805         jz=iz-1
     806         if (jz.lt.0) jz=ndz-1
    770807         jcl=jx+ndx*jy+nxy*jz+1
    771808         do ii=incell(jcl)+1,incell(jcl+1)
    … …  
    776813         jx=ix
    777814         jy=iy-1
     815         if (jy.lt.0) jy =ndy-1
    778816         jz=iz+1
     817         if (jz.ge.ndz) jz=0
    779818         jcl=jx+ndx*jy+nxy*jz+1
    780819         do ii=incell(jcl)+1,incell(jcl+1)
    … …  
    784823!        3
    785824         jx=ix-1
     825         if (jx.lt.0)jx =ndx-1
    786826         jy=iy+1
     827         if (jy.ge.ndy) jy=0
    787828         jz=iz
    788829         jcl=jx+ndx*jy+nxy*jz+1
    … …  
    793834!        4
    794835         jx=ix+1
     836         if (jx.ge.ndx) jx=0
    795837         jy=iy+1
     838         if (jy.ge.ndy) jy=0
    796839         jz=iz-1
     840         if (jz.lt.0) jz=ndz-1
    797841         jcl=jx+ndx*jy+nxy*jz+1
    798842         do ii=incell(jcl)+1,incell(jcl+1)
    … …  
    802846!        5
    803847         jx=ix+1
     848         if (jx.ge.ndx) jx = 0
    804849         jy=iy-1
     850         if (jy.lt.0) jy=ndy-1
    805851         jz=iz+1
     852         if (jz.ge.ndz) jz=0
    806853         jcl=jx+ndx*jy+nxy*jz+1
    807854         do ii=incell(jcl)+1,incell(jcl+1)
    … …  
    811858!        6
    812859         jx=ix+1
     860         if (jx.ge.ndx) jx=0
    813861         jy=iy-1
     862         if (jy.lt.0) jy=ndy-1
    814863         jz=iz-1
     864         if (jz.lt.0) jz=ndy-1
    815865         jcl=jx+ndx*jy+nxy*jz+1
    816866         do ii=incell(jcl)+1,incell(jcl+1)
    … …  
    895945     &                 +bsaexv(iatt1,iatt2)*r2
    896946                  etmp=etmp+etmp1
    897                   if (iat1.eq.43.and.iat2.eq.785) then
    898                      print *,'contribution ',iat1,iat2,etmp1
    899                      print *,'r2 = ',r2
    900                      print *,'atom types : ',iatt1,iatt2
    901                      print *,'sig2exv = ',sig2exv(iatt1,iatt2)
    902                      print *,'asa = ',asaexv(iatt1,iatt2)
    903                      print *,'bsa = ',bsaexv(iatt1,iatt2)
    904 
    905 !     call outpdb(1,'EXAMPLES/clash.pdb')
    906 !     stop
    907                   endif
    908947               else
    909948!                  print *,'atoms ', iat1,' and ',iat2,' were close',
Note: See TracChangeset for help on using the changeset viewer.