Changeset 32289cd


Ignore:
Timestamp:
11/19/09 11:29:41 (14 years ago)
Author:
baerbaer <baerbaer@…>
Branches:
master
Children:
38d77eb
Parents:
6650a56
Message:

Explicitly declare variables.

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

Files:
35 edited

Legend:

Unmodified
Added
Removed

  • INCL.H

    r6650a56 r32289cd  
    1 !      implicit none
    2       implicit integer*4 (i-n)
    3       implicit real*8 (a-h,o-z)
     1      implicit none
     2     
    43      character*255 version
    54

  • bgs.f

    r6650a56 r32289cd  
    22! Trial version implementing the semi-local conformational update
    33! BGS (Biased Gaussian Steps). This file presently contains the
    4 ! functions initlund, bgsposs and bgs. 
     4! functions initlund, bgsposs and bgs.
    55!
    66! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    … …  
    99
    1010
    11 ! Checks if it is possible to perform a BGS update starting at the 
     11! Checks if it is possible to perform a BGS update starting at the
    1212! variable indexed ipos. Calls: none.
    1313      logical function bgsposs(ips)
    1414      include 'INCL.H'
    1515      include 'incl_lund.h'
     16      integer jv, ips, iaa, nursvr, nnonfx, i
     17
    1618      logical ians
    1719
    … …  
    2022      ians=.true.
    2123!      print *,'evaluating bgs possibility for ',ips,nmvr(jv)
    22       if (nmvr(jv).ne.'phi') then 
     24      if (nmvr(jv).ne.'phi') then
    2325!         print *,'bgs not possible because variable name is ',nmvr(jv)
    2426         ians=.false.
    25       else if (iaa.gt.(irsml2(mlvr(jv))-3)) then 
     27      else if (iaa.gt.(irsml2(mlvr(jv))-3)) then
    2628!         print *,'bgs impossible, residue too close to end'
    2729!         print *,'iaa = ',iaa,' end = ',irsml2(mlvr(jv))
    2830         ians=.false.
    29       else 
    30          nnonfx=0 
     31      else
     32         nnonfx=0
    3133         do i=iaa,iaa+3
    32             if (iphi(i).gt.0) then 
     34            if (iphi(i).gt.0) then
    3335               if (.not.fxvr(iphi(i))) then
    3436                  nnonfx=nnonfx+1
    3537               endif
    3638            endif
    37             if (.not.fxvr(ipsi(i))) then 
     39            if (.not.fxvr(ipsi(i))) then
    3840               nnonfx=nnonfx+1
    3941            endif
    … …  
    5254
    5355! Biased Gaussian Steps. Implements a semi-local conformational update
    54 ! which modifies the protein backbone locally in a certain range of 
    55 ! amino acids. The 'down-stream' parts of the molecule outside the 
     56! which modifies the protein backbone locally in a certain range of
     57! amino acids. The 'down-stream' parts of the molecule outside the
    5658! region of update get small rigid body translations and rotations.
    5759!
    58 ! Use the update sparingly. It is rather local, and is not of great 
     60! Use the update sparingly. It is rather local, and is not of great
    5961! value if we need big changes in the conformation. It is recommended
    60 ! that this update be used to refine a structure around a low energy 
     62! that this update be used to refine a structure around a low energy
    6163! candidate structure. Even at low energies, if you always
    62 ! perform BGS, the chances of coming out of that minimum are small. 
    63 ! So, there is a probability bgsprob, which decides whether BGS or the 
    64 ! normal single angle update is used. 
     64! perform BGS, the chances of coming out of that minimum are small.
     65! So, there is a probability bgsprob, which decides whether BGS or the
     66! normal single angle update is used.
    6567!
    6668! Calls: energy, dummy (function provided as argument), addang, (rand)
    … …  
    7072      include 'INCL.H'
    7173      include 'incl_lund.h'
     74      double precision grnd, xiv, bv, ab, rv, dv, a, sum, p, r1, r2
     75      double precision ppsi, wfw, addang, enw, energy, wbw, rd, delta
     76      double precision dummy, eol1
     77
     78      integer ivar, i, nph, jv, ia, nursvr, icurraa, j, k, l
     79
    7280      external dummy
    7381      dimension xiv(8,3),bv(8,3),rv(3,3),dv(3,8,3)
    … …  
    7684! Initialize
    7785!      print *,'using BGS on angle ',nmvr(idvr(ivar))
    78       if (bgsnvar.eq.0) then 
     86      if (bgsnvar.eq.0) then
    7987         bgs=0
    8088         goto 171
    … …  
    9199      do i=1,4
    92100         icurraa=ia+i-1
    93          if (iphi(icurraa).gt.0.and..not.fxvr(iphi(icurraa))) then 
     101         if (iphi(icurraa).gt.0.and..not.fxvr(iphi(icurraa))) then
    94102            nph=nph+1
    95103            xiv(nph,1)=xat(iCa(icurraa))
    … …  
    102110     &           +bv(nph,3)*bv(nph,3)
    103111            iph(nph)=iphi(icurraa)
    104          endif 
     112         endif
    105113         if (.not.fxvr(ipsi(icurraa))) then
    106114            nph=nph+1
    … …  
    136144         enddo
    137145      enddo
    138       do i=1,nph 
     146      do i=1,nph
    139147         do j=i,nph
    140148            A(i,j)=0
    … …  
    157165               sum=sum-A(i,k)*A(j,k)
    158166            enddo
    159             if (i.eq.j) then 
    160                p(i)=sqrt(sum) 
    161             else 
     167            if (i.eq.j) then
     168               p(i)=sqrt(sum)
     169            else
    162170               A(j,i)=sum/p(i)
    163171            endif
    … …  
    177185         dph(i)=0
    178186      enddo
    179 ! Solve lower triangular matrix to get dphi proposals 
     187! Solve lower triangular matrix to get dphi proposals
    180188      do i=nph,1,-1
    181189         sum=ppsi(i)
    … …  
    190198      do i=1,nph
    191199         sum=sum+ppsi(i)*ppsi(i)
    192       enddo 
     200      enddo
    193201      wfw=exp(-sum)
    194202      do i=1,nph
    … …  
    208216      do i=1,4
    209217         icurraa=ia+i-1
    210          if (iphi(icurraa).gt.0.and..not.fxvr(iphi(icurraa))) then 
     218         if (iphi(icurraa).gt.0.and..not.fxvr(iphi(icurraa))) then
    211219            nph=nph+1
    212220            xiv(nph,1)=xat(iCa(icurraa))
    … …  
    219227     &           +bv(nph,3)*bv(nph,3)
    220228            iph(nph)=iphi(icurraa)
    221          endif 
    222          if (.not.fxvr(ipsi(icurraa))) then 
     229         endif
     230         if (.not.fxvr(ipsi(icurraa))) then
    223231            nph=nph+1
    224232            xiv(nph,1)=xat(iC(icurraa))
    … …  
    253261         enddo
    254262      enddo
    255       do i=1,nph 
     263      do i=1,nph
    256264         do j=i,nph
    257265            A(i,j)=0
    … …  
    262270            enddo
    263271            A(i,j)=bbgs*A(i,j)
    264             if (i.eq.j) then 
     272            if (i.eq.j) then
    265273               A(i,j)=A(i,j)+1
    266274            endif
    … …  
    274282               sum=sum-A(i,k)*A(j,k)
    275283            enddo
    276             if (i.eq.j) then 
    277                p(i)=sqrt(sum) 
    278             else 
     284            if (i.eq.j) then
     285               p(i)=sqrt(sum)
     286            else
    279287               A(j,i)=sum/p(i)
    280288            endif
    … …  
    285293         do j=i+1,nph
    286294            ppsi(i)=ppsi(i)+A(j,i)*dph(j)
    287          enddo 
     295         enddo
    288296      enddo
    289297      sum=0
    … …  
    291299         sum=sum+ppsi(i)*ppsi(i)
    292300      enddo
    293       wbw=exp(-sum) 
     301      wbw=exp(-sum)
    294302      do i=1,nph
    295303         wbw=wbw*p(i)
    … …  
    306314!       call outpdb(0,'after.pdb')
    307315!      print *,'after outpdb for after.pdb'
    308 !      do i=1,nph 
     316!      do i=1,nph
    309317!         print *,'BGS>',i,iph(i),vlvr(iph(i)),dph(i)
    310318!      enddo
    311       if (rd.ge.delta) then 
     319      if (rd.ge.delta) then
    312320!     accept
    313321         eol1=enw
    314322         bgs=1
    315323!         print *,'BGS move accepted'
    316       else 
     324      else
    317325!     reject
    318326         vlvr = ovr

  • cnteny.f

    r6650a56 r32289cd  
    44!
    55! Copyright 2005       Frank Eisenmenger, U.H.E. Hansmann,
    6 !                      Shura Hayryan, Chin-Ku 
     6!                      Shura Hayryan, Chin-Ku
    77! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    88!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    2020      include 'INCL.H'
    2121
    22       parameter (coeycn=2.d0)  ! min. cont. energy to display
     22      double precision coeycn, coey, eyatcn, cqi, xi, yi, zi, xij, yij
     23      double precision zij, rij2, rij4, rij6, rij, sr, ep, ey
     24
     25      integer ieyel, ntlvr, nml, iat1, nat, i, ifivr, i1s, io, iv, ia
     26      integer it, ic, i2s, ims, i1, i2, ii, ity, ivw, j, jj, jty, i14
     27      integer nbc, ir, nursat, idxat
     28
     29      parameter (coeycn=2.d0)  ! min. cont. energy to display
    2330
    2431      dimension eyatcn(mxat),idxat(mxat)

  • difang.f

    r6650a56 r32289cd  
    2121! ......................................................
    2222
    23       implicit real*8 (a-h,o-z)
    24 
     23      implicit none
     24      double precision pi, pi2, d, a2, a1
    2525      parameter (pi=3.141592653589793d0,
    2626     &           pi2=2.d0*pi)
    … …  
    4747! ......................................................
    4848
    49       implicit real*8 (a-h,o-z)
     49      implicit none
     50      double precision pi, pi2, d, a1, a2
    5051
    5152      parameter (pi=3.141592653589793d0,

  • energy.f

    r6650a56 r32289cd  
    44!
    55! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    6 !                      Shura Hayryan, Chin-Ku 
     6!                      Shura Hayryan, Chin-Ku
    77! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    88!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    2222
    2323      include 'INCL.H'
     24      double precision esm, teysl, enyshe, enyflx, enylun, enyreg
     25      double precision enysol, esolan, exvlun, eyabgn, eninteract
     26
     27      integer i
     28
    2429      double precision teysm, teyel, teyvw, teyhb, teyvr
    25 !      print *,'energy function with ientyp  = ',ientyp   
     30!      print *,'energy function with ientyp  = ',ientyp
    2631      esm = 0.d0
    2732      teysm = 0.d0
    … …  
    3237      teysl = 0.d0
    3338
    34       do i = 1,ntlml 
     39      do i = 1,ntlml
    3540         eysm=0
    3641         eyel=0
    … …  
    4449        if (ientyp.eq.0.or.ientyp.eq.3) then
    4550          esm=esm+enyshe(i)
    46         else if (ientyp.eq.1) then 
     51        else if (ientyp.eq.1) then
    4752          esm=esm+enyflx(i)
    4853        else if (ientyp.eq.2) then
    … …  
    5762!     The Lund term stores the hydrophobicity energy in eysl
    5863           teysl = teysl + eysl
    59         else 
     64        else
    6065!     .. and the excluded volume term in eyvw, which is calculated once.
    6166           teyvw = teyvw + eyvw
    6267        endif
    6368
    64         if (ireg.eq.1)  eyrg=enyreg(i)       
     69        if (ireg.eq.1)  eyrg=enyreg(i)
    6570
    6671      enddo
    6772
    68       if (ientyp.ne.2) then 
    69 !     Don't touch eysl if using Lund potential, as enylun stores 
     73      if (ientyp.ne.2) then
     74!     Don't touch eysl if using Lund potential, as enylun stores
    7075!     its hydrophobicity term there.
    7176         if (itysol.gt.0) then
    … …  
    7984            eysl=0.d0
    8085         endif
    81       else 
     86      else
    8287!     Add excluded volume term and save it in eyvw
    8388         esm=esm+exvlun(0)
    … …  
    8590      endif
    8691
    87 ! The Abagyan entropic corrections depend on the area exposed to the 
     92! The Abagyan entropic corrections depend on the area exposed to the
    8893! solvent for each residue. So, this term has to be evaluated after the
    8994! solvent term.
    9095      eyab=0.0
    9196      if (ientyp.eq.3) then
    92          do i = 1,ntlml 
     97         do i = 1,ntlml
    9398            eyab=eyab+eyabgn(i)
    9499         enddo
    … …  
    103108      eyvr = teyvr
    104109      eysl = teysl
    105      
     110
    106111      if (ientyp.ne.2) then
    107112!     This is temporary. eninteract() does not yet know how to calculate
    … …  
    116121!c Calculates the internal energy for a single molecule.
    117122!  All the partial energies are thus set to their values for molecule
    118 !  nml. 
     123!  nml.
    119124!
    120125!  @param nml the ID of the molecule
    … …  
    125130
    126131!f2py intent(in) nml
    127      
     132
    128133      include 'INCL.H'
     134      double precision esm, enyshe, enyflx, enylun, enyreg, enysol
     135      double precision esolan, exvlun, eyabgn
     136
     137      integer nml, i
     138
    129139      esm = 0.d0
    130140
    131       call setvar(nml,vlvr) 
     141      call setvar(nml,vlvr)
    132142
    133143      if (ientyp.eq.0.or.ientyp.eq.3) then
    … …  
    149159            eysl=0.d0
    150160         endif
    151       else 
     161      else
    152162         esm=esm+exvlun(nml)
    153163      endif
    154 ! The Abagyan entropic corrections depend on the area exposed to the 
     164! The Abagyan entropic corrections depend on the area exposed to the
    155165! solvent for each residue. So, this term has to be evaluated after the
    156166! solvent term.

  • eninteract.f

    r6650a56 r32289cd  
    1515
    1616        include 'INCL.H'
    17                
     17
     18      double precision cqi, xi, yi, zi, xj, yj, zj, xij, yij, zij, rij2
     19      double precision rij4, rij6, rij, sr, ep
     20
     21      integer iml, jml, ires, jres, iat, ity, jat, jty
     22
    1823        eysmi=0.0
    1924        eyeli=0.0
    … …  
    4247                    yj=yat(jat)
    4348                    zj=zat(jat)
    44                    
     49
    4550                    xij=xat(jat)-xi
    4651                    yij=yat(jat)-yi
    … …  
    6974     &                          -chb(ity,jty)/(rij6*rij4)
    7075                    endif
    71                   enddo 
     76                  enddo
    7277                enddo
    7378              enddo
    … …  
    7883        eysmi = eyeli + eyvwi + eyhbi
    7984        eninteract = eysmi
    80         return 
    81       end   
     85        return
     86      end

  • enylun.f

    r6650a56 r32289cd  
    11! *******************************************************************
    22! SMMP version of Anders Irback's force field, to be called the Lund
    3 ! force field. This file contains the function enylun, which in turn 
    4 ! calls all the terms in the energy function. The terms Bias (ebias), 
    5 ! Hydrogen bonds (ehbmm and ehbms), Hydrophobicity (ehp) and the 
    6 ! Excluded volume (eexvol and eloexv) are also implemented in this 
     3! force field. This file contains the function enylun, which in turn
     4! calls all the terms in the energy function. The terms Bias (ebias),
     5! Hydrogen bonds (ehbmm and ehbms), Hydrophobicity (ehp) and the
     6! Excluded volume (eexvol and eloexv) are also implemented in this
    77! file.
    88!
    … …  
    1313      include 'INCL.H'
    1414      include 'incl_lund.h'
     15      integer iat1, iat2
     16
    1517          integer ires,frst,lst,root
    1618          do iat1=iCa(ires)+frst,iCa(ires)+lst
    … …  
    3032      include 'INCL.H'
    3133      include 'incl_lund.h'
     34      double precision eunit, csacc
     35
     36      integer i, j, iml, iat1, iat2, k, iat3, l, iat4, m, iat3p, irs
     37      integer iatoff, iatmrg, ilp, lci
     38
    3239      character mynm*4
    3340      logical prlvr
    … …  
    8390            prlvr=.true.        ! residue i is a proline variant
    8491            print *, 'proline variant ',mynm,i
    85          else 
     92         else
    8693            prlvr=.false.
    87          endif 
    88 
    89          if (mynm.eq.'ala') then 
     94         endif
     95
     96         if (mynm.eq.'ala') then
    9097            nhpat(i)=1
    9198            ihpat(i,1)=iCa(i)+2
    … …  
    193200
    194201!     Initialization of the connections matrix matcon(i,j). The index
    195 !     i runs from -mxconr to +mxconr, and j from 1 to mxat. 
     202!     i runs from -mxconr to +mxconr, and j from 1 to mxat.
    196203!     matcon(i2-i1,i1) = 0, if the distance between atoms i1 and i2 is fixed
    197204!                      = 2, if atoms i1 and i2 are separated by 3 covalent
    198205!                           bonds and their distance can change
    199206!                      = 1, for all other pairs
    200 !     if abs(i2-i1) > mxconr, the atoms are assumed to be separated by 
    201 !     many bonds, and with no restriction on their distances. On a protein 
     207!     if abs(i2-i1) > mxconr, the atoms are assumed to be separated by
     208!     many bonds, and with no restriction on their distances. On a protein
    202209!     molecule made of natural amino acids, atoms with indices separated
    203210!     by more than 35 can not be connected by three covalent bonds.
    … …  
    228235                     do m=1,nbdat(iat2)
    229236                        iat3p=ibdat(m,iat2)
    230                         if (iat3p.ne.iat3) then 
     237                        if (iat3p.ne.iat3) then
    231238                           matcon(iat4-iat3p,iat3p)=2 ! 3 covalent bonds
    232239                           matcon(iat3p-iat4,iat4)=2 !
    … …  
    255262                matcon(iat2-iat1,iat1)=0
    256263                matcon(iat1-iat2,iat2)=0
    257             enddo 
     264            enddo
    258265            matcon(iat1-iCa(irsml1(iml))-1,iCa(irsml1(iml))+1)=2
    259266            matcon(iCa(irsml1(iml))+1-iat1,iat1)=2
    … …  
    266273!     Below: for certain residues, some atoms separated by 3 or more bonds
    267274!     do not change distance. So, the connection matrix term for such pairs
    268 !     should be zero. 
     275!     should be zero.
    269276
    270277         do irs=irsml1(iml),irsml2(iml)
    271278            if (irs.eq.irsml1(iml)) then
    272279               iatoff=1
    273             else 
     280            else
    274281               iatoff=0
    275282            endif
    … …  
    283290            if ((mynm.eq.'pro')) then
    284291               prlvr=.true.     ! residue i is a proline variant
    285             else 
     292            else
    286293               prlvr=.false.
    287             endif 
    288             if ((mynm.eq.'asn')) then 
     294            endif
     295            if ((mynm.eq.'asn')) then
    289296                call set_local_constr(irs,5,9,2)
    290             else if ((mynm.eq.'gln')) then 
     297            else if ((mynm.eq.'gln')) then
    291298                call set_local_constr(irs,8,12,5)
    292299            else if ((mynm.eq.'arg')) then
    … …  
    294301            else if ((mynm.eq.'his')) then
    295302                call set_local_constr(irs,5,12,2)
    296             else if (mynm.eq.'phe') then 
     303            else if (mynm.eq.'phe') then
    297304                call set_local_constr(irs,5,15,2)
    298305            else if (mynm.eq.'tyr') then
    … …  
    308315                call set_local_constr(irs,5,19,2)
    309316            else if (prlvr) then
    310 !           Proline. Many more distances are fixed because of the fixed 
     317!           Proline. Many more distances are fixed because of the fixed
    311318!           phi angle
    312319                call set_local_constr(irs,-3,11,-3)
    … …  
    347354         do iat1=iatrs1(irsml1(iml)),iatrs2(irsml2(iml))
    348355            do iat2=iat1+1,iatrs2(irsml2(iml))
    349                if (iat2-iat1.le.mxconr) then
    350                   if (matcon(iat2-iat1,iat1).eq.2) then
    351                     ilp=ilp+1
    352                     lcp1(ilp)=iat1
    353                     lcp2(ilp)=iat2
    354                   endif
     356               if ((iat2-iat1.le.mxconr).and.
     357     &                 matcon(iat2-iat1,iat1).eq.2) then
     358                  ilp=ilp+1
     359                  lcp1(ilp)=iat1
     360                  lcp2(ilp)=iat2
    355361               endif
    356362            enddo
    … …  
    358364
    359365         ilpnd(iml)=ilp
    360          if (iml.lt.ntlml) then 
     366         if (iml.lt.ntlml) then
    361367            ilpst(iml+1)=ilp+1
    362368         endif
    … …  
    369375         enddo
    370376      enddo
     377
    371378      print *,'finished initializing Lund force field'
    372379      end
    … …  
    374381      integer function ihptype(iaa)
    375382      include 'INCL.H'
     383      integer iaa, ityp
     384
    376385      character mynm*4
    377386      mynm=seq(iaa)
    … …  
    383392     &        .or.(mynm.eq.'met').or.(mynm.eq.'pro')) then
    384393         ityp=2
    385       else if ((mynm.eq.'phe').or.(mynm.eq.'tyr').or.(mynm.eq.'trp')) 
     394      else if ((mynm.eq.'phe').or.(mynm.eq.'tyr').or.(mynm.eq.'trp'))
    386395     &        then
    387396         ityp=3
    388397      endif
    389398      ihptype=ityp
    390       return 
     399      return
    391400      end
    392401
    … …  
    394403      include 'INCL.H'
    395404      include 'incl_lund.h'
     405      double precision q1, q2, et, xij, yij, zij
     406
     407      integer i, iat1, irsd, j, iat2
     408
    396409      dimension q1(2),q2(2)
    397410      data q1/-0.2,0.2/
    398       data q2/0.42,-0.42/ 
     411      data q2/0.42,-0.42/
    399412
    400413      et=0.0
    … …  
    412425      return
    413426      end
    414 !     Evaluates backbone backbone hydrogen bond strength for residues 
     427!     Evaluates backbone backbone hydrogen bond strength for residues
    415428!     i and j, taking the donor from residue i and acceptor from residue j
    416429      real*8 function ehbmmrs(i,j)
    417430      include 'INCL.H'
    418431      include 'incl_lund.h'
     432      double precision rdon2, racc2, evlu
     433
     434      integer i, j
     435
    419436      double precision r2,r4,r6,dx,dy,dz,ca,cb
    420437      integer d1,d2,a1,a2
    … …  
    427444      dz=zat(a1)-zat(d2)
    428445      r2=dx*dx+dy*dy+dz*dz
    429       if (r2.gt.cthb2) then 
     446      if (r2.gt.cthb2) then
    430447!         print *,'hbmm = 0 ',cthb2,r2,a1,a2,d1,d2
    431448!         print *,'a1,a2,d1,d2,r2 = ',a1,a2,d1,d2,r2,sighb2,cthb
    … …  
    464481      include 'INCL.H'
    465482      include 'incl_lund.h'
     483      double precision ebiasrs, shbm1, shbm2, etmp, ehbmmrs, ehp, etmp1
     484      double precision xij, yij, zij, r2, r6
     485
     486      integer istres, nml, indres, i, j, ihpi, ihptype, ihpj, i1, i2
     487      integer iat1, iat2, iatt1, iatt2
     488
    466489      character mynm*4
    467490      logical prlvr
    … …  
    485508     &           .or.(mynm.eq.'pron').or.(mynm.eq.'pro+')) then
    486509            prlvr=.true.        ! residue i is a proline variant
    487          else 
     510         else
    488511            prlvr=.false.
    489          endif 
     512         endif
    490513
    491514!     Bias, or local electrostatic term. Excluded from the list are
    … …  
    502525!     j run over the whole set of amino acids.
    503526!     No terms for residue i, if it is a proline variant.
    504          if (.not.prlvr) then 
     527         if (.not.prlvr) then
    505528            do j=istres,indres
    506                if ((j.lt.(i-2)).or.(j.gt.(i+1))) then 
     529               if ((j.lt.(i-2)).or.(j.gt.(i+1))) then
    507530                    shbm2=1.0
    508531                    if ((j.eq.istres).or.(j.eq.indres)) shbm2=0.5
    509532                    etmp=ehbmmrs(i,j)
    510533                    eyhb=eyhb+shbm1*shbm2*etmp
    511                 endif 
     534                endif
    512535            enddo
    513536         endif
    514537!     Hydrophobicity, only if residue i is hydrophobic to start with
    515538         ihpi=ihptype(i)
    516          if (ihpi.ge.0) then 
     539         if (ihpi.ge.0) then
    517540!        Unlike hydrogen bonds, the hydrophobicity potential is symmetric
    518541!        in i and j. So, the loop for j runs from i+1 to the end.
    … …  
    526549!                  print *, 'hydrophobicity contribution ',etmp,i,j
    527550                  eysl=eysl+etmp
    528                endif 
     551               endif
    529552            enddo
    530553         endif
    … …  
    535558!     Local pair or third-neighbour excluded volume
    536559!     Numerically this is normally the term with largest positive
    537 !     contribution to the energy in an equilibrated stystem. 
     560!     contribution to the energy in an equilibrated stystem.
    538561
    539562      i1=ilpst(nml)
    … …  
    551574         r2=(xij*xij + yij*yij + zij*zij)
    552575         if (r2.le.exvcut2) then
    553             r6=sig2lcp(iatt1,iatt2)/r2 
     576            r6=sig2lcp(iatt1,iatt2)/r2
    554577            r6=r6*r6*r6
    555578            etmp1=(r6*r6+asalcp(iatt1,iatt2)+bsalcp(iatt1,iatt2)*r2)
    … …  
    560583            etmp=etmp+etmp1
    561584         endif
    562 !         print *,'pair : ',iat1,iat2,' contribution ',etmp1 
     585!         print *,'pair : ',iat1,iat2,' contribution ',etmp1
    563586!         print *,exvcut2,r2
    564587      enddo
    … …  
    579602      include 'INCL.H'
    580603      include 'incl_lund.h'
     604      double precision b2, a2, r2min, xij, yij, zij, dtmp, ssum
     605
     606      integer ihp1, ihp2, ni, i1, nj, i2, i, k, j, l
     607
    581608      dimension r2min(12)
    582609
    … …  
    608635      ssum=0
    609636      do i=1,ni+nj
    610          if (r2min(i).le.b2) then 
    611             if (r2min(i).lt.a2) then 
     637         if (r2min(i).le.b2) then
     638            if (r2min(i).lt.a2) then
    612639               ssum=ssum+1
    613             else 
     640            else
    614641               ssum=ssum+(b2-r2min(i))/(b2-a2)
    615642            endif
    … …  
    626653      include 'incl_lund.h'
    627654!     For multi-chain systems it makes little sense to split the calculation
    628 !     of this term into an 'interaction part' and a contribution from 
     655!     of this term into an 'interaction part' and a contribution from
    629656!     individual molecules. So, normally this should always be called with
    630657!     argument nml=0. Only for diagnostic reasons, you might want to find
    631 !     the contribution from one molecule in a multi-chain system assuming 
     658!     the contribution from one molecule in a multi-chain system assuming
    632659!     there was no other molecule.
     660      double precision xmin, ymin, zmin, xmax, ymax, zmax, sizex, sizey
     661      double precision sizez, shiftx, shifty, shiftz, etmp, xij, yij
     662      double precision zij, r2, r6, etmp1
     663
     664      integer nml, istat, indat, i, incell, ndx, ndy, ndz, nxy, ncell
     665      integer nocccl, locccl, j, mx, my, mz, icellj, icell, jj, isort
     666      integer icl, lcell, ix, iz, iy, nex, ney, nez, nsame, nngbr, jx
     667      integer jy, jz, jcl, ii, ngbr, i1, iat1, i2, iat2, iatt1, iatt2
     668
    633669      dimension isort(mxat),ngbr(mxat),locccl(mxat),incell(mxcell)
    634670      dimension icell(mxat)
    635671      integer :: jx1, jy1, jz1
    636       logical doit
    637 
    638       if (nml.eq.0) then
     672
     673      if (nml.eq.0) then
    639674         istat=iatrs1(irsml1(1))
    640675         indat=iatrs2(irsml2(ntlml))
    641       else 
     676      else
    642677         istat=iatrs1(irsml1(nml))
    643678         indat=iatrs2(irsml2(nml))
    … …  
    645680
    646681      eyvw=0.0
    647 !     The beginning part of this implementation is very similar to the 
    648 !     assignment of cells to the atoms during calculation of solvent 
    649 !     accessible surface area. So, much of that part is similar. But 
     682!     The beginning part of this implementation is very similar to the
     683!     assignment of cells to the atoms during calculation of solvent
     684!     accessible surface area. So, much of that part is similar. But
    650685!     unlike the accessible surface calculations, this term is symmetric
    651686!     in any two participating atoms. So, the part after the assignment
    … …  
    724759         icellj=mx+my*ndx+mz*nxy+1
    725760         icell(j)=icellj
    726          if (icellj.gt.mxcell) then 
     761         if (icellj.gt.mxcell) then
    727762            print *,'exvlun> bad cell index ',icellj,' for atom ',j
    728763            stop
    729          else 
    730             if (incell(icellj).eq.0) then 
     764         else
     765            if (incell(icellj).eq.0) then
    731766!           previously unoccupied cell
    732767               nocccl=nocccl+1
    … …  
    756791         ix=mod(lcell-1,ndx)
    757792         iz = (lcell - 1) / nxy
    758          iy = ((lcell - 1) - iz * nxy) / ndx 
     793         iy = ((lcell - 1) - iz * nxy) / ndx
    759794
    760795c         print *,'icl=',icl,'absolute index of cell = ',lcell
    … …  
    767802         nngbr=0
    768803         do jx=ix,nex
    769             if (jx.ge.ndx) then 
     804            if (jx.ge.ndx) then
    770805               jx1=0
    771             else 
     806            else
    772807               jx1=jx
    773808            endif
    774809            do jy=iy,ney
    775                if (jy.ge.ndy) then 
     810               if (jy.ge.ndy) then
    776811                  jy1=0
    777                else 
     812               else
    778813                  jy1=jy
    779                endif 
     814               endif
    780815               do jz=iz,nez
    781                   if (jz.ge.ndz) then 
     816                  if (jz.ge.ndz) then
    782817                     jz1=0
    783                   else 
     818                  else
    784819                     jz1=jz
    785820                  endif
    786821                  jcl=jx1 + ndx*jy1 + nxy*jz1 + 1
    787 !                  write(*,*)'jcl,jx1,jy1,jz1:', jcl,jx1,jy1,jz1
     822!                  write (logString, *)'jcl,jx1,jy1,jz1:', jcl,jx1,jy1,jz1
    788823                  do ii=incell(jcl)+1,incell(jcl+1)
    789824!     count the total number of neighbours
    … …  
    799834         enddo
    800835!     A few more cells need to be searched, so that we cover 13 of the 26
    801 !     neighbouring cells. 
     836!     neighbouring cells.
    802837!        1
    803838         jx=ix+1
    804          if (jx.ge.ndx) jx=0 
     839         if (jx.ge.ndx) jx=0
    805840         jy=iy
    806841         jz=iz-1
    … …  
    883918               r2=(xij*xij+yij*yij+zij*zij)
    884919
    885                if (r2.le.exvcutg2) then
    886                   doit=.false.
    887                   if (abs(iat2-iat1).gt.mxconr ) then
    888                     doit=.true.
    889                   else if (matcon(iat2-iat1,iat1).eq.1) then
    890                     doit=.true.
    891                   endif
    892                   if (doit) then
     920               if (r2.le.exvcutg2) then
     921                  if (abs(iat2-iat1).gt.mxconr.or.
     922     &                 matcon(iat2-iat1,iat1).eq.1) then
    893923                     iatt1=ityat(iat1)
    894924                     iatt2=ityat(iat2)
    … …  
    931961      include 'incl_lund.h'
    932962
     963      double precision etmp, etmp1, xij, yij, zij, r2, r6
     964
     965      integer iat1, iat2, iatt1, iatt2
     966
    933967      etmp=0.0
    934968      etmp1=0.0
    … …  
    941975            r2=(xij*xij+yij*yij+zij*zij)
    942976
    943             if (r2.le.exvcutg2) then 
     977            if (r2.le.exvcutg2) then
    944978               if (abs(iat2-iat1).gt.mxconr.or.
    945979     &                 matcon(iat2-iat1,iat1).eq.1) then
    … …  
    951985     &                 +bsaexv(iatt1,iatt2)*r2
    952986                  etmp=etmp+etmp1
    953                else 
     987               else
    954988!                  print *,'atoms ', iat1,' and ',iat2,' were close',
    955989!     #                 'but matcon is ',matcon(iat2-iat1,iat1)

  • enyreg.f

    r6650a56 r32289cd  
    2121
    2222
     23      double precision eny
     24
     25      integer i, nml, j
     26
    2327      eny = 0.d0
    2428

  • enyshe.f

    r6650a56 r32289cd  
    11! **************************************************************
    22!
    3 ! This file contains the subroutines: enyshe 
     3! This file contains the subroutines: enyshe
    44!
    55! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    6 !                      Shura Hayryan, Chin-Ku 
     6!                      Shura Hayryan, Chin-Ku
    77! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    88!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    2020!
    2121! The function loops over all moving sets within the molecule. Within
    22 ! this loop it loops over the van-der-Waals domains of each atom in the 
     22! this loop it loops over the van-der-Waals domains of each atom in the
    2323! moving set and finally over the atoms that belong to the 1-4 interaction
    2424! set.
    … …  
    2828
    2929! If nml == 0 calculate the interaction between all pairs.
     30      double precision e0, vr, cqi, xi, yi, zi, xij, yij, zij, rij2
     31      double precision rij4, rij6, rij, sr, ep
     32
     33      integer nml, ntlvr, ifivr, i1s, io, iv, ia, it, ic, i2s, ims, i1
     34      integer i2, i, ity, ivw, j, jty, i14
     35
    3036      if (nml.eq.0) then
    3137          ntlvr = nvr
    … …  
    3339        ntlvr=nvrml(nml)
    3440      endif
    35      
     41
    3642      if (ntlvr.eq.0) then
    3743        write (*,'(a,i4)')
    … …  
    4854        ifivr = ivrml1(1)
    4955        i1s = imsml1(ntlml) + nmsml(ntlml)
    50       else 
     56      else
    5157! Index of first variable in molecule.
    5258        ifivr=ivrml1(nml)
    … …  
    5460        i1s=imsml1(nml)+nmsml(nml)
    5561      endif
    56 ! Loop over moving sets/variables in reverse order     
    57       do io=ifivr+ntlvr-1,ifivr,-1 
     62! Loop over moving sets/variables in reverse order
     63      do io=ifivr+ntlvr-1,ifivr,-1
    5864! The array iorvr contains the variables in an "apropriate" order.
    59         iv=iorvr(io)       
     65        iv=iorvr(io)
    6066! Index of the primary moving atom for the variable with index iv
    61         ia=iatvr(iv)       
     67        ia=iatvr(iv)
    6268! Get the type of variable iv (valence length, valence angle, dihedral angle)
    63         it=ityvr(iv)       
     69        it=ityvr(iv)
    6470! Class of variable iv's potential  (Q: What are they)
    65         ic=iclvr(iv)       
     71        ic=iclvr(iv)
    6672! If iv is a dihedral angle ...
    67         if (it.eq.3) then     
     73        if (it.eq.3) then
    6874! Barrier height * 1/2 of the potential of iv.
    6975          e0=e0to(ic)
    7076! Calculate the periodic potential term. sgto is the sign of the barrier, rnto is
    7177! the periodicity and toat is torsion angle(?) associate with atom ia.
    72           if (e0.ne.0.) 
     78          if (e0.ne.0.)
    7379     &         eyvr=eyvr+e0*(1.0+sgto(ic)*cos(toat(ia)*rnto(ic)))
    7480! else if iv is a valence angle ...
    75         elseif (it.eq.2) then 
     81        elseif (it.eq.2) then
    7682! vr is the valence angle of ia
    7783          vr=baat(ia)
    7884! else if iv is a valence length...
    79         elseif (it.eq.1) then 
     85        elseif (it.eq.1) then
    8086! vr is the length of the valence bond
    8187          vr=blat(ia)
    … …  
    8692        i2s=i1s-1
    8793! index of first moving set associated with iv
    88         i1s=imsvr1(iv) 
     94        i1s=imsvr1(iv)
    8995! Loop over all moving sets starting from the one associated with vr to the end.
    90         do ims=i1s,i2s 
     96        do ims=i1s,i2s
    9197! First atom of the current moving set
    9298          i1=latms1(ims)
    … …  
    94100          i2=latms2(ims)
    95101! Loop over all atoms of the current moving set.
    96           do i=i1,i2 
     102          do i=i1,i2
    97103! Atom class of current atom
    98104            ity=ityat(i)
    … …  
    104110            zi=zat(i)
    105111! Loop over the atoms of the van der Waals domain belonging to atom i
    106             do ivw=ivwat1(i),ivwat2(i) 
    107 ! Loop over the atoms of the van der Waals domain of the atoms of the 
     112            do ivw=ivwat1(i),ivwat2(i)
     113! Loop over the atoms of the van der Waals domain of the atoms of the
    108114! van der Waals domain of atom i
    109115! Q: Which atoms are in these domains?
    110               do j=lvwat1(ivw),lvwat2(ivw) 
     116              do j=lvwat1(ivw),lvwat2(ivw)
    111117! Atom type of partner
    112118                jty=ityat(j)
    … …  
    139145                endif
    140146
    141               enddo 
    142             enddo 
    143            
     147              enddo
     148            enddo
     149
    144150! Loop over 1-4 interaction partners
    145151! The interactions between atoms that are three bonds apart in the protein are
    146152! dominated by quantum mechanical effects. They are treated separately.
    147             do i14=i14at1(i),i14at2(i)   
     153            do i14=i14at1(i),i14at2(i)
    148154              j=l14at(i14)
    149155

  • enyshe_p.f

    r6650a56 r32289cd  
    11! **************************************************************
    22!
    3 ! This file contains the subroutines: enyshe 
     3! This file contains the subroutines: enyshe
    44!
    55! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    6 !                      Shura Hayryan, Chin-Ku 
     6!                      Shura Hayryan, Chin-Ku
    77! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    88!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    1414
    1515!     ............................................................................
    16 !     
     16!
    1717!     PURPOSE: Calculate internal energy of molecule 'nml' with ECEPP parameters
    18 !     
     18!
    1919!     CALLS: none
    20 !     
     20!
    2121!     The function loops over all moving sets within the molecule. Within
    22 !     this loop it loops over the van-der-Waals domains of each atom in the 
     22!     this loop it loops over the van-der-Waals domains of each atom in the
    2323!     moving set and finally over the atoms that belong to the 1-4 interaction
    2424!     set.
    … …  
    3131
    3232! If nml == 0 calculate the interaction between all pairs.
     33      double precision eysmsum, teysm, teyel, teyvw, teyhb, teyvr
     34      double precision startwtime, e0, vr, cqi, xi, yi, zi, xij, yij
     35      double precision zij, rij2, rij4, rij6, rij, sr, ep, endwtime
     36
     37      integer nml, ntlvr, ifivr, i1s, iend, istart, loopcounter, io, iv
     38      integer ia, it, ic, i2s, ims, i1, i2, i, ity, ivw, j, jty, i14
     39      integer ierror
     40
    3341      if (nml.eq.0) then
    3442         ntlvr = nvr
    … …  
    3644         ntlvr=nvrml(nml)
    3745      endif
    38      
     46
    3947      if (ntlvr.eq.0) then
    4048         write (*,'(a,i4)')
    … …  
    5664         ifivr = ivrml1(1)
    5765         i1s = imsml1(ntlml) + nmsml(ntlml)
    58       else 
     66      else
    5967!     Index of first variable in molecule.
    6068         ifivr=ivrml1(nml)
    … …  
    6270         i1s=imsml1(nml)+nmsml(nml)
    6371      endif
    64 !     Loop over variables in reverse order     
     72!     Loop over variables in reverse order
    6573!     This is the first loop to parallize. We'll just split the moving sets
    6674!     over the number of available processors and sum the energy up in the end.
    … …  
    7280      startwtime = MPI_Wtime()
    7381      loopcounter = 0
    74 !      do io=ifivr+ntlvr-1,ifivr,-1 
    75       do io = workPerProcessor(nml, myrank) - 1, 
     82!      do io=ifivr+ntlvr-1,ifivr,-1
     83      do io = workPerProcessor(nml, myrank) - 1,
    7684     &        workPerProcessor(nml, myrank+1), -1
    7785         if (io.lt.istart) then
    … …  
    7987         endif
    8088!     The array iorvr contains the variables in an "apropriate" order.
    81          iv=iorvr(io)       
     89         iv=iorvr(io)
    8290!     Index of the primary moving atom for the variable with index iv
    83          ia=iatvr(iv)       
     91         ia=iatvr(iv)
    8492!     Get the type of variable iv (valence length, valence angle, dihedral angle)
    85          it=ityvr(iv)       
     93         it=ityvr(iv)
    8694!     Class of variable iv's potential  (Q: What are they)
    87          ic=iclvr(iv)       
     95         ic=iclvr(iv)
    8896!     If iv is a dihedral angle ...
    89          if (it.eq.3) then     
     97         if (it.eq.3) then
    9098!     Barrier height * 1/2 of the potential of iv.
    9199            e0=e0to(ic)
    92100!     Calculate the periodic potential term. sgto is the sign of the barrier, rnto is
    93101!     the periodicity and toat is torsion angle(?) associate with atom ia.
    94             if (e0.ne.0.) 
     102            if (e0.ne.0.)
    95103     &           teyvr=teyvr+e0*(1.0+sgto(ic)*cos(toat(ia)*rnto(ic)))
    96104!     else if iv is a valence angle ...
    97          elseif (it.eq.2) then 
     105         elseif (it.eq.2) then
    98106!     vr is the valence angle of ia
    99107            vr=baat(ia)
    100108!     else if iv is a valence length...
    101          elseif (it.eq.1) then 
     109         elseif (it.eq.1) then
    102110!     vr is the length of the valence bond
    103111            vr=blat(ia)
    … …  
    108116         i2s=i1s-1
    109117!     index of first moving set associated with iv
    110          i1s=imsvr1(iv) 
     118         i1s=imsvr1(iv)
    111119!     Loop over all moving sets starting from the one associated with vr to the end.
    112          do ims=i1s,i2s 
     120         do ims=i1s,i2s
    113121!     First atom of the current moving set
    114122            i1=latms1(ims)
    … …  
    116124            i2=latms2(ims)
    117125!     Loop over all atoms of the current moving set.
    118             do i=i1,i2 
     126            do i=i1,i2
    119127!     Atom class of current atom
    120128               ity=ityat(i)
    … …  
    126134               zi=zat(i)
    127135!     Loop over the atoms of the van der Waals domain belonging to atom i
    128                do ivw=ivwat1(i),ivwat2(i) 
    129 !     Loop over the atoms of the van der Waals domain of the atoms of the 
     136               do ivw=ivwat1(i),ivwat2(i)
     137!     Loop over the atoms of the van der Waals domain of the atoms of the
    130138!     van der Waals domain of atom i
    131139!     Q: Which atoms are in these domains?
    132                   do j=lvwat1(ivw),lvwat2(ivw) 
     140                  do j=lvwat1(ivw),lvwat2(ivw)
    133141
    134142                     loopcounter = loopcounter + 1
    … …  
    163171                     endif
    164172
    165                   enddo 
    166                enddo 
    167                
     173                  enddo
     174               enddo
     175
    168176!     Loop over 1-4 interaction partners
    169177!     The interactions between atoms that are three bonds apart in the protein are
    170178!     dominated by quantum mechanical effects. They are treated separately.
    171                do i14=i14at1(i),i14at2(i)   
     179               do i14=i14at1(i),i14at2(i)
    172180                  loopcounter = loopcounter + 1
    173181                  j=l14at(i14)
    … …  
    210218
    211219!     Collect energies from all nodes and sum them up
    212       call MPI_ALLREDUCE(teysm, eysmsum, 1, MPI_DOUBLE_PRECISION, 
     220      call MPI_ALLREDUCE(teysm, eysmsum, 1, MPI_DOUBLE_PRECISION,
    213221     &     MPI_SUM, my_mpi_comm, ierror)
    214222      call MPI_ALLREDUCE(teyel, eyel, 1, MPI_DOUBLE_PRECISION, MPI_SUM,
    215223     &     my_mpi_comm, ierror)
    216       call MPI_ALLREDUCE(teyvw, eyvw, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 
    217      &     my_mpi_comm, ierror)
    218       call MPI_ALLREDUCE(teyhb, eyhb, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 
    219      &     my_mpi_comm, ierror)
    220       call MPI_ALLREDUCE(teyvr, eyvr, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 
     224      call MPI_ALLREDUCE(teyvw, eyvw, 1, MPI_DOUBLE_PRECISION, MPI_SUM,
     225     &     my_mpi_comm, ierror)
     226      call MPI_ALLREDUCE(teyhb, eyhb, 1, MPI_DOUBLE_PRECISION, MPI_SUM,
     227     &     my_mpi_comm, ierror)
     228      call MPI_ALLREDUCE(teyvr, eyvr, 1, MPI_DOUBLE_PRECISION, MPI_SUM,
    221229     &     my_mpi_comm, ierror)
    222230

  • enysol.f

    r6650a56 r32289cd  
    44!
    55! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    6 !                      Shura Hayryan, Chin-Ku 
     6!                      Shura Hayryan, Chin-Ku
    77! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    88!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    1010! **************************************************************
    1111
    12      
     12
    1313      real*8 function enysol(nmol)
    1414
    … …  
    1717!
    1818!     Double Cubic Lattice algorithm for calculating the
    19 !     solvation energy of proteins using 
     19!     solvation energy of proteins using
    2020!     solvent accessible area method.
    2121!
    … …  
    2424!
    2525! -------------------------------------------------------------
    26 ! TODO: Check the solvent energy for multiple molecules     
     26! TODO: Check the solvent energy for multiple molecules
    2727!     arguments
    2828      integer nmol
    29      
     29
    3030!     functions
    3131      integer nursat
    … …  
    4141      double precision sizes, trad, zmin, xmax, xmin, ymax, zmax
    4242      double precision sdr, sdd, volume
    43      
     43
    4444      dimension numbox(mxat),inbox(mxbox+1),indsort(mxat),look(mxat)
    4545      dimension xyz(mxinbox,3),radb(mxinbox),radb2(mxinbox)
    … …  
    4747
    4848!      common/ressurf/surfres(mxrs)
    49      
     49
    5050      eyslh = 0.0
    5151      eyslp = 0.0
    … …  
    6161      do i=nrslow,nrshi
    6262       surfres(i) = 0.0d0
    63       end do 
     63      end do
    6464
    6565      numat= nup - nlow + 1
    … …  
    6868         inbox(i)=0
    6969      end do
    70      
     70
    7171      asa=0.0d0
    7272      vdvol=0.0d0
    … …  
    135135       stop
    136136      end if
    137        
     137
    138138! Let us shift the borders to home all boxes
    139139
    … …  
    175175        inbox(i+1)=inbox(i+1)+inbox(i)
    176176      end do
    177          
    178        
     177
     178
    179179!   Sorting the atoms by the their box numbers
    180180
    … …  
    184184         indsort(jj)=i
    185185         inbox(j)=jj-1
    186       end do   
    187          
     186      end do
     187
    188188! Getting started
    189189
    … …  
    205205             ney=min(iy+1,ndy-1)
    206206             nez=min(iz+1,ndz-1)
    207                      
     207
    208208!  Atoms in the boxes around
    209209
    … …  
    221221               end do
    222222              end do
    223              end do     
    224                              
     223             end do
     224
    225225             do  ia=inbox(ibox)+1,inbox(ibox+1)
    226226               jbi=indsort(ia)
    … …  
    274274                     end do
    275275 99                  continue
    276    
     276
    277277                     if(ik.gt.nnei)then
    278278                       surfc(il)=.true.
    … …  
    316316                 surfres(jres) = surfres(jres) + area
    317317               end if
    318              end do 
     318             end do
    319319           end if
    320320          end do
    … …  
    341341      character lin*80
    342342      integer i
    343 !    Skipping comment lines, which begin with '!' 
     343!    Skipping comment lines, which begin with '!'
    344344      read(20,'(a)') lin
    345345      do while(lin(1:1).eq.'!')
    … …  
    352352!       write(*,'(a,i5)') 'the number of points---->',npnt
    353353
    354 !    Read the surface points   
     354!    Read the surface points
    355355
    356356      do i=1,npnt
    … …  
    358358!          write(31,'(3f20.10)') spoint(i,1),spoint(i,2),spoint(i,3)
    359359      end do
    360  
     360
    361361      return
    362  
     362
    363363      end
    364364

  • enysol_p.f

    r6650a56 r32289cd  
    44!
    55! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    6 !                      Shura Hayryan, Chin-Ku 
     6!                      Shura Hayryan, Chin-Ku
    77! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    88!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    1010! **************************************************************
    1111
    12      
     12
    1313      real*8 function enysol(nmol)
    1414
    … …  
    1919!
    2020!     Double Cubic Lattice algorithm for calculating the
    21 !     solvation energy of proteins using 
     21!     solvation energy of proteins using
    2222!     solvent accessible area method.
    2323!
    … …  
    2626!
    2727! -------------------------------------------------------------
    28 ! TODO: Check the solvent energy for multiple molecules     
     28! TODO: Check the solvent energy for multiple molecules
     29      double precision startwtime, avr_x, avr_y, avr_z, xmin, ymin, zmin
     30      double precision xmax, ymax, zmax, rmax, diamax, sizex, sizey
     31      double precision sizez, shiftx, shifty, shiftz, boxpp, trad, dx
     32      double precision dy, dz, dd, akrad, dr, xyz, radb, radb2, sdd, sdr
     33      double precision area, volume, eyslsum, endwtime
     34
     35      integer nmol, nrslow, nrshi, nlow, nup, i, numat, inbox, j, ndx
     36      integer ndy, ndz, nqxy, ncbox, mx, my, mz, nboxj, numbox, jj
     37      integer indsort, iboxmin, iboxmax, ibox, iz, iy, ix, lbn, nsx, nsy
     38      integer nsz, nex, ney, nez, jcnt, jz, jy, jx, jbox, ii, look, ia
     39      integer jbi, nnei, ib, jtk, il, lst, ilk, ik, icount, jres, nursat
     40      integer ierror, nhx, mhx, nbt, mbt
     41
    2942      dimension numbox(mxat),inbox(mxbox+1),indsort(mxat),look(mxat)
    3043      dimension xyz(mxinbox,3),radb(mxinbox),radb2(mxinbox)
    … …  
    3346
    3447!       common/ressurf/surfres(mxrs)
    35       real*8 tsurfres(mxrs) 
     48      real*8 tsurfres(mxrs)
    3649      startwtime = MPI_Wtime()
    3750      root = 0
    … …  
    5063      do i=nrslow,nrshi
    5164       surfres(i) = 0.0d0
    52       end do 
     65      end do
    5366
    5467      numat= nup - nlow + 1
    … …  
    5770         inbox(i)=0
    5871      end do
    59      
     72
    6073      asa=0.0d0
    6174      vdvol=0.0d0
    … …  
    124137       stop
    125138      end if
    126        
     139
    127140! Let us shift the borders to home all boxes
    128141
    … …  
    164177        inbox(i+1)=inbox(i+1)+inbox(i)
    165178      end do
    166          
    167        
     179
     180
    168181!   Sorting the atoms by the their box numbers
    169182
    … …  
    173186         indsort(jj)=i
    174187         inbox(j)=jj-1
    175       end do   
    176          
     188      end do
     189
    177190!    Getting started
    178191!    We have to loop over ncbox boxes and have no processors available
    … …  
    202215             ney=min(iy+1,ndy-1)
    203216             nez=min(iz+1,ndz-1)
    204                      
     217
    205218!  Atoms in the boxes around
    206219
    … …  
    218231               end do
    219232              end do
    220              end do     
    221                              
     233             end do
     234
    222235             do  ia=inbox(ibox)+1,inbox(ibox+1)
    223236               jbi=indsort(ia)
    … …  
    271284                     end do
    272285 99                  continue
    273    
     286
    274287                     if(ik.gt.nnei)then
    275288                       surfc(il)=.true.
    … …  
    313326                 surfres(jres) = surfres(jres) + area
    314327               end if
    315              end do 
     328             end do
    316329           end if
    317330!           end do
    318331!          end do
    319332       end do
    320       call MPI_ALLREDUCE(eysl, eyslsum, 1, MPI_DOUBLE_PRECISION, 
     333      call MPI_ALLREDUCE(eysl, eyslsum, 1, MPI_DOUBLE_PRECISION,
    321334     &      MPI_SUM,my_mpi_comm, ierror)
    322335!       write(*,*) 'enysol>', myrank, eysl, eyslsum
    323336      tsurfres = surfres
    324       call MPI_ALLREDUCE(tsurfres, surfres, mxrs, MPI_DOUBLE_PRECISION, 
     337      call MPI_ALLREDUCE(tsurfres, surfres, mxrs, MPI_DOUBLE_PRECISION,
    325338     &      MPI_SUM,my_mpi_comm, ierror)
    326339      eysl = eyslsum
    327      
     340
    328341      endwtime = MPI_Wtime()
    329342      if (myrank.le.-1) then
    … …  
    349362      subroutine tessel
    350363      include 'INCL.H'
     364      integer i
     365
    351366      character lin*80
    352367
    353 !    Skipping comment lines, which begin with '!' 
     368!    Skipping comment lines, which begin with '!'
    354369
    355370      read(20,'(a)') lin
    … …  
    363378!        write(*,'(a,i5)') 'the number of points---->',npnt
    364379
    365 !    Read the surface points   
     380!    Read the surface points
    366381
    367382      do i=1,npnt
    368383         read(20,'(3f20.10)') spoint(i,1),spoint(i,2),spoint(i,3)
    369          
     384
    370385!        write(31,'(3f20.10)') spoint(i,1),spoint(i,2),spoint(i,3)
    371386      end do
    372  
     387
    373388      return
    374  
     389
    375390      end
    376391

  • esolan.f

    r6650a56 r32289cd  
    44!
    55! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    6 !                      Shura Hayryan, Chin-Ku 
     6!                      Shura Hayryan, Chin-Ku
    77! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    88!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    1414! -----------------------------------------------------------------
    1515!
    16 ! Calculates the solvation energy of the protein with 
    17 ! solavtion parameters model E=\sum\sigma_iA_i. 
     16! Calculates the solvation energy of the protein with
     17! solavtion parameters model E=\sum\sigma_iA_i.
    1818! The solvent accessible surface area per atom
    19 ! and its gradients with respect to the Cartesian coordinates of 
     19! and its gradients with respect to the Cartesian coordinates of
    2020! the atoms are calculated using the projection method described in
    2121!
    … …  
    2525!
    2626! INPUT: nmol - the order number of the protein chain.
    27 !      The atomic coordinates, radii and solvation parameters 
     27!      The atomic coordinates, radii and solvation parameters
    2828!      are taken from the global arrays of SMMP
    2929!
    30 ! OUTPUT: global array gradan(mxat,3) contains the gradients of 
     30! OUTPUT: global array gradan(mxat,3) contains the gradients of
    3131!         solvation energy per atom
    3232!         local array as(mxat) contains accessible surface area per atom
    … …  
    4242
    4343      include 'INCL.H'
    44      
    45      
     44
     45
    4646      integer nmol, ks0, ks2
    4747Cf2py intent(in) nmol
    48      
     48
    4949      parameter (ks0=mxat*(mxat+1))
    5050      parameter (ks2=mxat+mxat)
    … …  
    5353
    5454      integer neib, neibp, ivrx, iv, ial, i, ij, j, iii, idi,ivr, ifu,
    55      &          ita2, ita3, ine, kk, j2, jkk, jjj, jj, jk, k, l, ll, 
     55     &          ita2, ita3, ine, kk, j2, jkk, jjj, jj, jk, k, l, ll,
    5656     &          nb, numat, nyx
    5757      real*8 vertex, ax, pol, as, ayx, ayx1, probe, dd, ddat, dadx, gp,
    5858     &          dalp,dbet, daalp, dabet, vrx, dv, dx, dy, dz, dt, di,
    59      &          dii, ss, dta, dtb, di1, di2, gs, xold, yold, zold, 
     59     &          dii, ss, dta, dtb, di1, di2, gs, xold, yold, zold,
    6060     &          energy, sss, a1, a, ab1, a2, aa, ac2, ac1, ab2, am2,
    6161     &          am, aia, ak, alp, am1, ap1, amibe, amabe, amaal, amial,
    … …  
    6363     &          ba, bcd, bc, bet, c1, c, c2, cc, sfcg, caba, cg, cfsg,
    6464     &          cfcg, daba, ct, cs, d1, d2, dbe, dal, dcotbma, dcbet,
    65      &          dcalp, dcbetmalp, dcbetpalp, ddp2, div, dsbetmalp, 
    66      &          dsalp, dsbet, yy, dsbetpalp, enr, p1, p2, r42, r22, 
     65     &          dcalp, dcbetmalp, dcbetpalp, ddp2, div, dsbetmalp,
     66     &          dsalp, dsbet, yy, dsbetpalp, enr, p1, p2, r42, r22,
    6767     &          r22p, r2, pom, prob, r, r42p, ratom, rr, s, sfsg, sf,
    68      &          sg, vv1, t, ufi, tt, uga, uv, vv, vv2, wv, xx, xv, yv, 
     68     &          sg, vv1, t, ufi, tt, uga, uv, vv, vv2, wv, xx, xv, yv,
    6969     &          zz, zv
    7070      dimension neib(0:mxat),neibp(0:mxat),ivrx(ks0)
    … …  
    7878      integer ta2(0:mxat),ta3(0:mxat),fullarc(0:mxat),al(0:ks2)
    7979      real*8 neibor(mxat,4)
    80      
     80
    8181      real*8, dimension(:, :, :), allocatable :: grad
    82      
     82
    8383      allocate(grad(mxat, mxat, 3))
    84      
    85      
     84
     85
    8686!     open(3,file='solvation.dat')! Output file with solvation data
    8787
    … …  
    9090      if (nmol.eq.0) then
    9191         numat=iatrs2(irsml2(ntlml))-iatrs1(irsml1(1))+1 ! Number of atoms
    92       else 
     92      else
    9393         numat=iatrs2(irsml2(nmol))-iatrs1(irsml1(nmol))+1 ! Number of atoms
    9494      endif
    … …  
    9696      do i=1,numat
    9797        do j=1,3
    98          gradan(i,j)=0.0d0 
     98         gradan(i,j)=0.0d0
    9999      end do
    100100      end do
    … …  
    102102      do i=1,numat
    103103         xold(i)=xat(i)! Redefine the coordinates just for safety
    104          yold(i)=yat(i)! 
    105          zold(i)=zat(i)! 
    106       pol(i)=rvdw(i)*rvdw(i)!The water radius is already added in 'main' 
     104         yold(i)=yat(i)!
     105         zold(i)=zat(i)!
     106      pol(i)=rvdw(i)*rvdw(i)!The water radius is already added in 'main'
    107107      As(i)=pi4*pol(i)! Initially the whole surface of the atom is accessible.
    108108      end do
    109      
     109
    110110             do iii=1,numat
    111111              do jjj=1,numat
    … …  
    120120!***************************
    121121
    122       do 1 i=1,numat ! Lop over atoms "i"           
     122      do 1 i=1,numat ! Lop over atoms "i"
    123123! ----------------------------------------------------
    124          
    125        R=rvdw(i) 
     124
     125       R=rvdw(i)
    126126       jj=0
    127127        do j=1,numat !Find the neighbours of "i"
    … …  
    132132         write(*,*)'ERROR in data: centres of two atoms coincide!'
    133133         write(*,*)i,j,xold(i),yold(i),zold(i),rvdw(i),
    134      &                   xold(j),yold(j),zold(j),rvdw(j) 
     134     &                   xold(j),yold(j),zold(j),rvdw(j)
    135135         stop 'Centres of atoms coincide!'
    136136       endif
    … …  
    146146         end if
    147147         end do!Neighbours
    148        
     148
    149149       neib(0)=jj   ! The number of neighbors of "i"
    150150
    … …  
    197197!
    198198        uga=grnd() ! Random \gamma angle
    199         ufi=grnd() ! Random \pi angle 
     199        ufi=grnd() ! Random \pi angle
    200200        uga=pi*uga/2
    201201        ufi=pi*ufi/2
    … …  
    219219        end do
    220220      end do
    221 !       
    222      
     221!
     222
    223223       if(jjj.ne.0) then ! Rotation is necessary
    224224         sfsg=sf*sg
    … …  
    243243      do jj=1,neib(0)
    244244      j=neib(jj)
    245        neibor(j,1)=(ddat(j,2))**2+(ddat(j,3))**2+     
     245       neibor(j,1)=(ddat(j,2))**2+(ddat(j,3))**2+
    246246     &            (ddat(j,4)-R)**2-pol(j)           ! a
    247247        neibor(j,2)=-pom*ddat(j,2)                  ! b
    … …  
    259259       do while(k.gt.1)               ! B
    260260        k=k-1
    261 ! Analyse mutual disposition of every pair of neighbours                                                     
     261! Analyse mutual disposition of every pair of neighbours
    262262        do 13 L=neib(0),k+1,-1        ! A
    263263
    … …  
    283283          elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).le.      ! a01 02
    284284     &    -dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
    285 ! The circle neib(k) encloses circle neib(L) and the later is discarded 
     285! The circle neib(k) encloses circle neib(L) and the later is discarded
    286286           neib(0)=neib(0)-1
    287287           do j=L,neib(0)
    … …  
    326326          elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).ge.  ! a03 02
    327327     &      dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
    328 ! Don't exclude neib(k) 
     328! Don't exclude neib(k)
    329329           neib(0)=neib(0)-1
    330330           do j=k,neib(0)
    … …  
    370370          elseif(D.le.0d0.and.dsqrt((b1-b2)**2+(c1-c2)**2).ge.   ! a07 02
    371371     &      dsqrt(b2*b2+c2*c2-4d0*d2)+dsqrt(b1*b1+c1*c1-4d0*d1)) then
    372 ! discard the circle neib(L) 
     372! discard the circle neib(L)
    373373           neib(0)=neib(0)-1
    374374           do j=L,neib(0)
    … …  
    385385           p2=(b1-b2)*pom
    386386!-- Assign t and s coordinates and the order number of circles
    387            vertex(iv,1)=(t+p1)/am 
     387           vertex(iv,1)=(t+p1)/am
    388388           vertex(iv,2)=(s-p2)/am
    389389           vertex(iv,3)=neib(k)
    … …  
    394394           vertex(iv,3)=neib(k)
    395395           vertex(iv,4)=neib(L)
    396 !-- 
     396!--
    397397          endif                                                 ! 10
    398398
    … …  
    491491            dbe=amabe-amibe
    492492            if(dal.lt.pi) then
    493              As(i)=R22*(pi-dal) 
     493             As(i)=R22*(pi-dal)
    494494            elseif(dbe.lt.pi) then
    495495             As(i)=R22*(pi-dbe)
    … …  
    510510       enddo
    511511       do j=1,neib(0)
    512 ! "iv" is the number of intersection point. 
     512! "iv" is the number of intersection point.
    513513        do k=1,iv
    514514       if(neibp(j).eq.vertex(k,3).or.neibp(j).eq.vertex(k,4)) goto 30
    … …  
    532532       do k=1,ifu
    533533         a=neibor(fullarc(k),1)
    534          if(a.lt.0.d0) then 
     534         if(a.lt.0.d0) then
    535535          aia=-1.d0
    536536         else
    … …  
    579579         ak=a*a
    580580         daba=dabs(a)
    581          if(a.lt.0d0) then 
     581         if(a.lt.0d0) then
    582582          aia=-1d0
    583583          else
    … …  
    679679          do j=1,nyx
    680680
    681 ! 
     681!
    682682!  Escape 'bad' intersections.
    683683           if(dabs(ayx1(j+1)-ayx1(j)).lt.1d-8) goto 40
    … …  
    686686           ap1=ct+rr*dcos(prob) ! The middle point of the arc.
    687687           ap2=cs+rr*dsin(prob)
    688 ! Verify if the middle point belongs to a covered part. 
     688! Verify if the middle point belongs to a covered part.
    689689! If yes then omit.
    690690           do ij=1,ial
    … …  
    754754             dabet(3)=(a*ac2-ak*c2+wv*daba*a2*dsbet)/am2
    755755             dabet(4)=(2d0*ak*a2)/am2
    756                dv(2)=R42*b*vv1     
    757                dv(3)=R42*c*vv1     
    758                dv(4)=(d-R42*a)*vv1     
    759                dv(1)=-dv(4)*R42     
    760                dx(1)=R42*vv1-(d+R42*a)*dv(1)*vv2     
    761              dx(2)=-(d+R42*a)*dv(2)*vv2     
    762              dx(3)=-(d+R42*a)*dv(3)*vv2     
    763              dx(4)=1d0*vv1-(d+R42*a)*dv(4)*vv2     
     756               dv(2)=R42*b*vv1
     757               dv(3)=R42*c*vv1
     758               dv(4)=(d-R42*a)*vv1
     759               dv(1)=-dv(4)*R42
     760               dx(1)=R42*vv1-(d+R42*a)*dv(1)*vv2
     761             dx(2)=-(d+R42*a)*dv(2)*vv2
     762             dx(3)=-(d+R42*a)*dv(3)*vv2
     763             dx(4)=1d0*vv1-(d+R42*a)*dv(4)*vv2
    764764             dy(1)=(-2d0*d+4d0*R42*a)/daba*vv1-(b*b+c*c-2d0*a*d+2d0*
    765765     &             R42*ak)*(aia*vv+daba*dv(1))/ak*vv2
    … …  
    769769     &             R42*ak)*aia*dv(3)/a*vv2
    770770               dy(4)=-2.d0*a/daba*vv1-(b*b+c*c-2.d0*a*d+2.d0*
    771      &             R42*ak)*aia*dv(4)/a*vv2     
     771     &             R42*ak)*aia*dv(4)/a*vv2
    772772             dz(1)=-2d0*d/wv/a*vv1-wv*(vv+a*dv(1))/ak*vv2
    773773             dz(2)=b/wv/a*vv1-wv*dv(2)/a*vv2
    … …  
    820820             dii(4,2)=dii(1,2)*R42
    821821             dii(4,3)=2d0*(ddat(jk,4)+R)*R42
    822              
     822
    823823             do idi=1,3
    824824              grad(i,jk,idi)=grad(i,jk,idi)+
    … …  
    826826     &                   di(3)*dii(3,idi)+di(4)*dii(4,idi)
    827827
    828              enddo             
     828             enddo
    829829             do idi=1,4
    830830               dta(idi)=5d-1*daalp(idi)*(((yv-zv*((-b*dsbetpalp
    … …  
    838838             di1(idi)=R2*(aia*daalp(idi)-4d0*xv*dta(idi)/(4d0+tt*tt))
    839839            di2(idi)=R2*(-aia*dabet(idi)-4d0*xv*dtb(idi)/(4d0+tt*tt))
    840              enddo             
     840             enddo
    841841
    842842             dii(1,1)=2d0*ddat(kk,2)
    … …  
    86686640           continue
    867867          enddo
    868          
     868
    869869         else ! analyse the case with lines (not necessary after rotation).
    870870
    … …  
    920920           enddo
    921921           probe(nyx+1)=ayx1(nyx)+1d0
    922            
     922
    923923
    924924           bc=b*b+c*c
    … …  
    96096022           continue
    961961            endif
    962          enddo       
     962         enddo
    96396311      continue
    964964
    … …  
    10171017111   continue
    10181018
    1019 !     write(3,*)'   No   Area    sigma   Enrg    gradx   grady',     
     1019!     write(3,*)'   No   Area    sigma   Enrg    gradx   grady',
    10201020!    & '   gradz   Rad    Atom'
    10211021!     write(3,*)
    10221022
    1023      
     1023
    10241024      do i=1,numat
    10251025        enr=As(i)*sigma(i)
    … …  
    10381038
    10391039      esolan=sss
    1040      
     1040
    10411041      deallocate(grad)
    10421042!     close(3)
    … …  
    10461046200   format(2i4,3f16.6)
    10471047700   format(i5,7f8.3,5x,a4)
    1048          
     1048
    10491049      end

  • eyabgn.f

    r6650a56 r32289cd  
    55!                      Jan H. Meinke, Sandipan Mohanty
    66!
    7 ! Corrections to ECEPP energy terms due to R. A. Abagyan et al. 
    8 ! 
     7! Corrections to ECEPP energy terms due to R. A. Abagyan et al.
     8!
    99! Two terms are calculated: eyrccr and eyentr, representing respectively
    10 ! c a term to slightly shift the backbone dihedral angle preferences in 
     10! c a term to slightly shift the backbone dihedral angle preferences in
    1111! the ECEPP potential slightly away from the helix region, and another
    12 ! term to estimate the side-chain entropy from a given configuration. 
     12! term to estimate the side-chain entropy from a given configuration.
    1313!
    1414!
    … …  
    1616      real*8 function eyrccr(nml)
    1717      include 'INCL.H'
     18      double precision et, rsscl
     19
     20      integer nml, istres, indres, i, ipsi, in, ica, ic, mlvr, iphi
     21
    1822      dimension iN(mxrs),iCa(mxrs),iC(mxrs),mlvr(mxvr)
    1923      dimension iphi(mxrs),ipsi(mxrs)
    … …  
    3438         call tolost(mynm)
    3539         if ((mynm.eq.'val').or.(mynm.eq.'ile').or.
    36      &           (mynm.eq.'thr')) then 
     40     &           (mynm.eq.'thr')) then
    3741            rsscl=1.0
    38          else 
     42         else
    3943            rsscl=0.5
    4044         endif
    … …  
    4852      eyrccr=et
    4953      return
    50       end 
     54      end
    5155
    5256      subroutine init_abgn
    … …  
    5458!       dimension rsstrg(mxrs)
    5559!       common /abgncor/rsstrg
     60      double precision xarea, estrg
     61
     62      integer i, istres, indres, imytyp, j
     63
    5664      dimension xarea(nrsty),estrg(nrsty)
    5765      character mynm*4
    … …  
    102110!            print *,'comparing ',mynm,' with ',rsnmcd(j),imytyp
    103111         enddo
    104          if (imytyp.eq.0) then 
     112         if (imytyp.eq.0) then
    105113            print *, 'Unknown residue type ',seq(i)
    106114            print *, 'Abagyan term strength set to 0'
    107115            rsstrg(i)=0.0
    108          else 
     116         else
    109117            rsstrg(i)=estrg(imytyp)/xarea(imytyp)
    110118         endif
    111119!         print *,'residue ',i,seq(i),' type ',imytyp
    112 !         print *, 'strength for residue ',i,seq(i),' is ',rsstrg(i) 
     120!         print *, 'strength for residue ',i,seq(i),' is ',rsstrg(i)
    113121      enddo
    114122      print *, 'initialized Abagyan corrections to ECEPP force field'
    115123      end
    116      
     124
    117125      real*8 function eyentr(nml)
    118126      include 'INCL.H'
    … …  
    121129!       common/ressurf/surfres(mxrs)
    122130!      common/bet/beta
     131      double precision eentr, aars, strh
     132
     133      integer nml, istres, indres, i
     134
    123135      eentr=0
    124136      if (nml.eq.0) then
    … …  
    136148!        The maximal burial entropies were estimated at temperature 300k
    137149!        The values in the array estrg are k_B * T (=300k) * Entropy
    138 !        Presently we need it at temperature 1/beta, so we need to 
     150!        Presently we need it at temperature 1/beta, so we need to
    139151!        multiply the strengths in estrg with (1/beta)/(300 kelvin)
    140152!        300 kelvin is approximately 0.59576607 kcal/mol.
    … …  
    149161      return
    150162      end
    151      
    152       real*8 function eyabgn(nml) 
     163
     164      real*8 function eyabgn(nml)
    153165      include 'INCL.H'
     166      double precision eyrccr, eyentr
     167
     168      integer nml
     169
    154170      eyabgn=eyrccr(nml)+eyentr(nml)
    155171!      print *,'Abagyan term = ',eyabgn

  • gradient.f

    r6650a56 r32289cd  
    44!
    55! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    6 !                      Shura Hayryan, Chin-Ku 
     6!                      Shura Hayryan, Chin-Ku
    77! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    88!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    2121      include 'INCL.H'
    2222
     23
     24      double precision esm
     25
     26      integer i, ivr1, ivr2, j
    2327
    2428      esm = 0.d0

  • init_energy.f

    r6650a56 r32289cd  
    55!
    66! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    7 !                      Shura Hayryan, Chin-Ku 
     7!                      Shura Hayryan, Chin-Ku
    88! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    99!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    1717! PURPOSE: initialize energy parameters
    1818!        0  => ECEPP2 or ECEPP3 depending on the value of sh2
    19 !        1  => FLEX 
     19!        1  => FLEX
    2020!        2  => Lund force field
    2121!        3  => ECEPP with Abagyan corrections
    … …  
    2929      include 'INCL.H'
    3030
     31      integer ll, iendst, its
     32
    3133      character libdir*(*),tesfil*80
    32      
     34
    3335      if (ientyp.eq.1) then
    3436          flex = .true.
    35       else 
     37      else
    3638          flex = .false.
    3739      end if
    38      
     40
    3941      lunlib=10
    4042      ll=iendst(libdir)
    … …  
    5658      endif
    5759
    58       if (ientyp .eq. 2) then 
     60      if (ientyp .eq. 2) then
    5961        reslib=libdir(1:ll)//'lib.lun'
    6062      endif
    … …  
    7880
    7981      else
    80  
     82
    8183        if (itysol.ne.0) then
    8284          write(*,'(a)') ' init_energy>  undefined solvent type !'
    … …  
    110112
    111113      include 'INCL.H'
     114
     115      double precision hbc, hba, ri, ai, aei, aic, a, aj, c, rij
     116
     117      integer i, j, iac, ido, jac, jdo
    112118
    113119      dimension hbc(mxhbdo,mxhbac),hba(mxhbdo,mxhbac)
    … …  
    308314
    309315      include 'INCL.H'
     316     
     317      integer i
    310318
    311319!  Atom types ------------------------------------------------------------
    … …  
    517525     & 'ile ','leu ','lys ','lys+','met ','phe ','cpro','pro ','cpru',
    518526     & 'prou','pron','pro+','ser ','thr ','trp ','tyr ','val ' /
    519  
     527
    520528      data (onltcd(i),i=1,nrsty)/  ! One-letter codes for amino acid types
    521529     & 'A',   'R',   'R',   'N',   'D',   'D',   'C',   'C',   'Q',
    … …  
    527535!      coefficients for their solvation free energy (kcal/molxA**2)
    528536
    529 !  Method: 
    530 
    531 !  itysol=1 : OONS --> T.Ooi, et al, 
     537!  Method:
     538
     539!  itysol=1 : OONS --> T.Ooi, et al,
    532540!                      Proc. Natl. Acad. Sci. USA 8 (1987) 3086-3090.
    533 !  itysol=2 : JRF  --> J.Vila, et al, 
     541!  itysol=2 : JRF  --> J.Vila, et al,
    534542!                      PROTEINS: Struct Funct Genet 10(1991) 199-218.
    535 !  itysol=3 : WE92 --> L.Wesson, D.Eisenberg, 
     543!  itysol=3 : WE92 --> L.Wesson, D.Eisenberg,
    536544!                      Protein Science 1 (1992) 227-235.
    537545!  itysol=4 : SCH1 --> D.Eisenberg, et al,
    … …  
    539547!  itysol=5 : SCH2 --> A.H.Juffer, et al,
    540548!                      Proteine Science 4 (1995) 2499-2509.
    541 !  itysol=6 : SCH3 --> L.Wesson, D.Eisenberg, 
     549!  itysol=6 : SCH3 --> L.Wesson, D.Eisenberg,
    542550!                      Protein Science 1 (1992) 227-235.
    543551!  itysol=7 : SCH4 --> C.A. Schiffer, et al,
    544552!                      Mol. Simul. 10(1993) 121-149.
    545 !  itysol=8 : EM86 --> D.Eisenberg, A.D. Mclachlan, 
     553!  itysol=8 : EM86 --> D.Eisenberg, A.D. Mclachlan,
    546554!                      Nature 319 (1986) 199-203.
    547555!  itysol=9 : BM   --> B. Freyberg, et al,

  • init_molecule.f

    r6650a56 r32289cd  
    44!
    55! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    6 !                      Shura Hayryan, Chin-Ku 
     6!                      Shura Hayryan, Chin-Ku
    77! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    88!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    1616! PURPOSE: construct starting structure of molecule(s)
    1717!
    18 !          iabin = 1  : ab Initio using sequence & 
     18!          iabin = 1  : ab Initio using sequence &
    1919!                       variables given in input files
    2020!          iabin != 1 : sequence, variable information
    … …  
    3232      include 'INCP.H'
    3333
     34      integer iabin, iendst, ntl, i, j, l, it, ier, ir, nursvr, i1, i2
     35      integer its
     36
    3437Cf2py character*80 optional, intent(in) :: seqfile = ' '
    35 Cf2py character*80 optional, intent(in) :: varfile = ' ' 
    36      
     38Cf2py character*80 optional, intent(in) :: varfile = ' '
     39
    3740      character grpn*4,grpc*4
    38       character navr*3, nars*4 
     41      character navr*3, nars*4
    3942      character seqfile*80, varfile*80
    4043      integer ontlml
    … …  
    4447      readFromStdin = .false.
    4548
    46       write (*,*) 'init_molecule: Solvent: ', itysol
    47       if (iabin.eq.1) then 
    48          
     49      write (logString, *) 'init_molecule: Solvent: ', itysol
     50      if (iabin.eq.1) then
     51
    4952!     ----------------------------------------- get sequence for molecule(s)
    5053         lunseq=11
    … …  
    5356         endif
    5457         if (iendst(seqfile).le.1.or.seqfile.eq.' ') then
    55  1          write (*,'(/,a,$)') ' file with SEQUENCE:'
     58 1          write (logString, '(/,a,$)') ' file with SEQUENCE:'
    5659            seqfil=' '
    5760            read (*,'(a)',err=1) seqfil
    5861            readFromStdin = .true.
    59          else 
     62         else
    6063            seqfil = seqfile
    61          endif 
     64         endif
    6265         call redseq
    63          
    64          write (*,*) 'File with sequence is ', seqfil(1:iendst(seqfil))
    65          
     66
     67         write (logString, *) 'File with sequence is ',
     68     &      seqfil(1:iendst(seqfil))
     69
    6670!     --------------------------------- read & assemble data from libraries
    6771!     initial coordinates, interaction lists
    68          
     72
    6973         ntl = ntlml
    7074         do i=ontlml, ntl
    71            
     75
    7276            call getmol(i)      ! assemble data from libraries
    73            
     77
    7478            do j=1,6            ! initialize global parameters
    7579               gbpr(j,i)=0.d0
    7680            enddo
    77            
     81
    7882            call bldmol(i)      ! co-ordinates
    79            
     83
    8084            ntlml = i
    8185            call addend(i,grpn,grpc) ! modify ends
    8286            call setmvs(i) ! determine sets of moving atoms for given variables
    8387            call mklist(i)      ! compile lists of interaction partners
    84            
     88
    8589         enddo
    86          
     90
    8791!     --------------------------- Read the initial conformation if necessary
    88          if(readFromStdin) then 
    89             write (*,'(a,$)') ' file with VARIABLES:'
    90 !     
     92         if(readFromStdin) then
     93            write (logString, '(a,$)') ' file with VARIABLES:'
     94!
    9195            varfil=' '
    9296            read(*,'(a)',end=2,err=2) varfil
    … …  
    96100         l=iendst(varfil)
    97101         if (l.gt.0.and.varfil.ne.' ') then
    98             write (*,'(1x,a,/)') varfil(1:l)
     102            write (logString, '(1x,a,/)') varfil(1:l)
    99103            lunvar=13
    100            
     104
    101105            call redvar         ! get vars. and rebuild
    102            
    103          endif
    104          
    105  2       write(*,*) ' '
    106          
    107          ireg = 0
    108          
    109       else                      ! =========================== from PDB
    110          if (iendst(seqfile).le.1) then
    111  3          write (*,'(/,a,$)') ' PDB-file:'
    112             seqfil=' '
    113             read (*,'(a)',err=3) seqfil
    114          else
    115             seqfil = seqfile
    116          endif
    117          write (*,*) 'PDB structure ',seqfil(1:iendst(seqfil))
    118          print *, 'calling readpdb with ',seqfile
    119          call pdbread(seqfil,ier)
    120          
    121          if (ier.ne.0) stop
    122          
    123          call pdbvars()
    124          
    125          ireg = 1
    126          
    127       endif
    128      
    129 ! If Lund force field is in use, keep omega angles fixed
    130       if (ientyp.eq.2) then
    131          do iv=1,nvrml(ntlml)
    132             if ((nmvr(iv)(1:2).eq.'om')) then
    133                 vlvr(iv)=pi
    134                 toat(iatvr(iv))=pi
    135                 fxvr(iv)=.true.
    136                 print *, 'Fixed variable ',iv,nmvr(iv),vlvr(iv)
    137             endif
    138          enddo
    139       endif
     106
     107         endif
     108
     109 2       write (logString, *) ' '
    140110
    141111!     -------------------- get: nvr,idvr, vlvr, olvlvr
    … …  
    161131         enddo
    162132
     133         ireg = 0
     134
     135      else                      ! =========================== from PDB
     136         if (iendst(seqfile).le.1) then
     137 3          write (logString, '(/,a,$)') ' PDB-file:'
     138            seqfil=' '
     139            read (*,'(a)',err=3) seqfil
     140         else
     141            seqfil = seqfile
     142         endif
     143         write (logString, *) 'PDB structure ',seqfil(1:iendst(seqfil))
     144         print *, 'calling readpdb with ',seqfile
     145         call pdbread(seqfil,ier)
     146
     147         if (ier.ne.0) stop
     148
     149         call pdbvars()
     150
     151         ireg = 1
     152
     153      endif
     154
    163155!     -------------------------- set var. amplitudes for simulations
    164      
     156
    165157      do i=1,ivrml1(ntlml)+nvrml(ntlml)-1
    166          
     158
    167159         if (ityvr(i).eq.3.and..not.fxvr(i)) then ! torsion
    168            
     160
    169161            navr = nmvr(i)
    170            
     162
    171163            ir = nursvr(i)
    172164            nars = seq(ir)
    173            
     165
    174166            if (                         navr(1:2).eq.'om'
    175            
     167
    176168     &     .or.nars(1:3).eq.'arg'.and.(navr(1:2).eq.'x5'
    177169     &           .or.navr(1:2).eq.'x6')
    178            
     170
    179171     &           .or.(nars(1:3).eq.'asn'.or.nars(1:3).eq.'asp')
    180172     &           .and.navr(1:2).eq.'x3'
    181            
     173
    182174     &           .or.(nars(1:3).eq.'gln'.or.nars(1:3).eq.'glu')
    183175     &           .and.navr(1:2).eq.'x4'
    184            
     176
    185177     &           ) then
    186            
     178
    187179!     axvr(i) = pi/9.d0  ! 20 deg.
    188180            axvr(i) = pi2       ! Trying out 360 deg. for these as well
    189            
     181
    190182         else
    191183            axvr(i) = pi2       ! 360 deg.
    192184         endif
    193          
     185
    194186      else
    195187         axvr(i) = 0.d0
    196188      endif
    197      
     189
    198190      enddo                     ! vars.
    199      
     191
    200192!     --------------------- initialize solvation pars. if necessary
    201193
    202194      if (itysol.ne.0) then
    203          
     195
    204196         i1=iatrs1(irsml1(1))   ! 1st atom of 1st molecule
    205197         i2=iatrs2(irsml2(ntlml)) ! last atom of last molecule
    206          
     198
    207199         its = iabs(itysol)
    208          
     200
    209201         do i=i1,i2             ! all atoms
    210202            it=ityat(i)
    211203            sigma(i)=coef_sl(its,it)
    212204            rvdw(i) =rad_vdw(its,it)
    213            
     205
    214206            if (nmat(i)(1:1).ne.'h') rvdw(i)=rvdw(i)+rwater
    215            
     207
    216208         enddo
    217          
     209
    218210      endif
    219211! Initialize calpha array
    220212      do i=ontlml, ntlml
    221          call c_alfa(i,1)   
     213         call c_alfa(i,1)
    222214      enddo
    223215
    224216!     Initialize arrays used in the BGS update
    225       call init_lund()     
     217      call init_lund()
    226218      return
    227219      end
    228      
    229      
     220
     221

  • main.f

    r6650a56 r32289cd  
    44!
    55! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    6 !                      Shura Hayryan, Chin-Ku 
     6!                      Shura Hayryan, Chin-Ku
    77! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    88!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    1616! $Id: main.f 334 2007-08-07 09:23:59Z meinke $
    1717! **************************************************************
    18      
     18
    1919      program main
    2020
    2121      include 'INCL.H'
    2222      include 'INCP.H'
     23      double precision eps, temp, e, energy
     24
     25      integer maxloglevel, logfileunit, iabin, imin, maxit, nequi
     26      integer nsweep, nmes, ncalls, nacalls
     27
    2328      common/updstats/ncalls(5),nacalls(5)
    2429      character*80 libdir, seqfile, varfile
    … …  
    3237!            Directory for SMMP libraries
    3338!     Change the following directory path to where you want to put SMMP
    34 !     libraries of residues. 
     39!     libraries of residues.
    3540      libdir='./SMMP/'
    36      
     41
    3742!     Set the maximum log level. The larger the number the more detailed
    3843!     the log.
    … …  
    4853      ientyp = 0
    4954!        0  => ECEPP2 or ECEPP3 depending on the value of sh2
    50 !        1  => FLEX 
     55!        1  => FLEX
    5156!        2  => Lund force field
    5257!        3  => ECEPP with Abagyan corrections
    … …  
    7075      iabin = 1  ! =0: read from PDB-file
    7176                 ! =1: ab Initio from sequence (& variables)
    72       seqfile='EXAMPLES/enkefa.seq'
    73       varfile='EXAMPLES/enkefa.ann'
    74 !       varfile = ' '
    75      
     77      seqfile='polyq.seq'
     78!      seqfile='polyA.pdb'
     79      varfile='polyq.var'
     80       varfile = ' '
     81
    7682      ntlml = 0
    7783      write (*,*) 'Solvent: ', itysol
    7884!     Initialize random number generator.
    7985      call sgrnd(31433)
    80      
     86
    8187      if (itysol.eq.0.and.ientyp.eq.3) then
    8288         print *,'Can not use Abagyan entropic corrections without '
    … …  
    8793      call init_molecule(iabin,grpn,grpc,seqfile,varfile)
    8894
    89 ! Decide if and when to use BGS, and initialize Lund data structures 
     95! Decide if and when to use BGS, and initialize Lund data structures
    9096      bgsprob=0.75   ! Prob for BGS, given that it is possible
    91 ! upchswitch= 0 => No BGS 1 => BGS with probability bgsprob 
    92 ! 2 => temperature dependent choice 
     97! upchswitch= 0 => No BGS 1 => BGS with probability bgsprob
     98! 2 => temperature dependent choice
    9399      upchswitch=1
    94100      rndord=.true.
    … …  
    96102      if (ientyp.eq.2) call init_lundff
    97103      if (ientyp.eq.3) call init_abgn
    98      
     104
    99105
    100106! ========================================  Add your task down here
    … …  
    103109      maxit = 15000 ! maximum number of iterations in minimization
    104110      eps = 1.0d-7 ! requested precision
    105       call minim(imin, maxit, eps)
    106       call outvar(0, ' ')
     111!      call minim(imin, maxit, eps)
     112!      call outvar(0, ' ')
    107113!     To do a canonical Monte Carlo simulation uncomment the lines below
    108 !       nequi = 100
    109 !       nsweep = 50000
    110 !       nmes = 10
    111 !       temp = 300.0
    112 !       lrand = .true.
     114       nequi = 100
     115       nsweep = 50000
     116       nmes = 10
     117       temp = 300.0
     118       lrand = .true.
     119       E = energy()
     120       write(*,*) E, eyel,eyvw,eyhb,eyvr
     121       call outpdb(1, "polyA.pdb")
    113122!      Canonical Monte Carlo
    114 !       call canon(nequi, nsweep, nmes, temp, lrand)
     123!        call canon(nequi, nsweep, nmes, temp, lrand)
    115124
    116125!      For simulated annealing uncomment the lines below
    … …  
    118127!      tmax = 500.0
    119128!      call anneal(nequi, nsweep, nmes, tmax, tmin, lrand);
    120 ! ========================================  End of main     
     129! ========================================  End of main
    121130       end

  • main_p.f

    r6650a56 r32289cd  
    11!     **************************************************************
    2 !     
     2!
    33!     This file contains the   main (PARALLEL TEMPERING  JOBS ONLY,
    44!     FOR SINGULAR PROCESSOR JOBS USE main)
    5 !     
     5!
    66!     This file contains also the subroutine: p_init_molecule
    7 !     
     7!
    88!     Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    9 !     Shura Hayryan, Chin-Ku 
     9!     Shura Hayryan, Chin-Ku
    1010! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    1111!                      Jan H. Meinke, Sandipan Mohanty
    12 !     
     12!
    1313!     CALLS init_energy,p_init_molecule,partem_p
    14 !     
     14!
    1515!     **************************************************************
    1616      program pmain
    … …  
    2121      include 'mpif.h'
    2222
     23      double precision startwtime, group_world, error, endwtime
     24
     25      integer ierr, num_proc, iabin, nequi, nswp, nmes, nsave, ifrm, j
     26      integer i, nml, nresi, my_pt_rank, ncalls, nacalls
     27
    2328      character*80 libdir
    2429      character*80 in_fil,ou_fil,filebase, varfile
    … …  
    2833      logical newsta
    2934
    30 !c    Number of replicas 
     35!c    Number of replicas
    3136      integer num_replica
    3237!c    Number of processors per replica
    … …  
    6267
    6368!     =================================================== Energy setup
    64       libdir='SMMP/'     
     69      libdir='SMMP/'
    6570!     Directory for SMMP libraries
    6671
    … …  
    7176      ientyp = 0
    7277!     0  => ECEPP2 or ECEPP3 depending on the value of sh2
    73 !     1  => FLEX 
     78!     1  => FLEX
    7479!     2  => Lund force field
    7580!     3  => ECEPP with Abagyan corrections
    76 !     
     81!
    7782
    7883      sh2=.false.               ! .true. for ECEPP/2; .false. for ECEPP3
    … …  
    106111                            !   temperatures must have at least as many
    107112                            !   entries
    108       nequi=10              ! Number of MC sweeps before measurements 
    109                             !   and replica exchanges are started 
     113      nequi=10              ! Number of MC sweeps before measurements
     114                            !   and replica exchanges are started
    110115      nswp=500000           ! Number of sweeps
    111116      nmes=10               ! Interval for measurements and replica exchange
    112117      nsave=1000            ! Not used at the moment
    113                    
    114       switch = -1           ! How should the configuration be   
     118
     119      switch = -1           ! How should the configuration be
    115120                            !   initialized?
    116                             ! -1 stretched chain 
     121                            ! -1 stretched chain
    117122                            !  0 don't do anything
    118123                            !  1 initialize each angle to a random value
    119      
     124
    120125      ifrm=0
    121126      ntlml = 0
    122127
    123 ! Decide if and when to use BGS, and initialize Lund data structures 
     128! Decide if and when to use BGS, and initialize Lund data structures
    124129      bgsprob=0.6    ! Prob for BGS, given that it is possible
    125 ! upchswitch= 0 => No BGS 1 => BGS with probability bgsprob 
    126 ! 2 => temperature dependent choice 
     130! upchswitch= 0 => No BGS 1 => BGS with probability bgsprob
     131! 2 => temperature dependent choice
    127132      upchswitch=1
    128133      rndord=.true.
    129134!     =================================================================
    130135!     Distribute nodes to parallel tempering tasks
    131 !     I assume that the number of nodes available is an integer 
     136!     I assume that the number of nodes available is an integer
    132137!     multiple n of the number of replicas. Each replica then gets n
    133138!     processors to do its energy calculation.
    … …  
    139144!     could just use i * num_ppr but this way it's more flexible.
    140145      j = 0
    141       do i = 1, num_replica 
    142          ranks(i) = j 
     146      do i = 1, num_replica
     147         ranks(i) = j
    143148         proc_range(1) = j
    144149         proc_range(2) = j + num_ppr - 1
    … …  
    146151         call mpi_group_range_incl(group_world, 1, proc_range, group(i)
    147152     &                              ,error)
    148          write (*,*) "Assigning rank ", j, proc_range, 
     153         write (*,*) "Assigning rank ", j, proc_range,
    149154     &               "to group", group(i)
    150155         call flush(6)
    … …  
    167172      call mpi_group_incl(group_world, num_replica, ranks, group_partem,
    168173     &                    error)
    169       call mpi_comm_create(mpi_comm_world, group_partem, partem_comm, 
     174      call mpi_comm_create(mpi_comm_world, group_partem, partem_comm,
    170175     &                     error)
    171176
    … …  
    184189         varfile = 'EXAMPLES/1bdd.var'
    185190         call init_molecule(iabin, grpn, grpc,in_fil,varfile)
    186       else 
     191      else
    187192         filebase = "conf_0000.var"
    188193         call init_molecule(iabin, grpn, grpc,in_fil,
    … …  
    220225      nml = 1
    221226      call distributeWorkLoad(no, nml)
    222      
     227
    223228      call partem_p(num_replica, nequi, nswp, nmes, nsave, newsta,
    224229     &              switch, rep_id, partem_comm)

  • metropolis.f

    r6650a56 r32289cd  
    2020      include 'INCP.H'
    2121      include 'incl_lund.h'
     22     
     23      double precision eol, energy, grnd, dv, addang, enw, delta, dummy
     24      double precision rd, ex, acz, eol1
     25
     26      integer nsw, iupstate, iupt, ivar, jv, iml, i, ncalls, nacalls
     27
    2228      external dummy
    2329!      common/bet/beta
    … …  
    178184     
    179185      subroutine accanalyze(iuptype,iupdstate)
     186      integer ncalls, nacalls, iupdstate, iuptype
    180187      common/updstats/ncalls(5),nacalls(5)
    181188      ncalls(5)=ncalls(5)+1
    … …  
    202209      integer iiii
    203210      double precision bbbb,curprob, grnd,up2bmax,up2bmin
     211
    204212      common/updtparam/up2bmax,up2bmin
    205213      curprob=(bbbb-up2bmin)/(up2bmax-up2bmin)
    … …  
    211219      end
    212220      block data updtchs
     221      integer ncalls, nacalls
    213222      double precision up2bmax, up2bmin
    214223      common/updtparam/up2bmax,up2bmin

  • mincjg.f

    r6650a56 r32289cd  
    44!
    55! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    6 !                      Shura Hayryan, Chin-Ku 
     6!                      Shura Hayryan, Chin-Ku
    77! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    88!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    4040! ....................................................................
    4141
    42       implicit real*8 (a-h,o-z)
    43       implicit integer*4 (i-n)
     42      double precision amf, tol, eps, gsqrd, fmin, f, t, g, ga, d, xa, x
     43      double precision gsq2, gnew, dfpr, stmin, fuit, gdit, gt, gdmi
     44      double precision sbound, stepch, step, wo, sum, fch, acur, ddspln
     45      double precision gspln, beta, gama, yt, dt, gamden, di
     46
     47      integer mxfcon, maxlin, mxn, ier, iter, iterfm, iterrs, nfun
     48      integer nfopt, i, n, nfbeg, iretry, maxfun
    4449
    4550      parameter (AMF = 10.d0,
    … …  
    102107
    103108    2 step = stmin + stepch
    104       wo = 0.d0 
     109      wo = 0.d0
    105110
    106111      do i=1,n
    … …  
    182187
    183188      if ( fch .ne. 0.d0 )  ddspln = ddspln + (wo + wo) / stepch
    184  
     189
    185190      if ( gdmi .eq. 0.d0 )  goto 6
    186191

  • nursvr.f

    r6650a56 r32289cd  
    44!
    55! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    6 !                      Shura Hayryan, Chin-Ku 
     6!                      Shura Hayryan, Chin-Ku
    77! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    88!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    1818! ...........................................................
    1919      include 'INCL.H'
     20
     21      integer i, ifirs, ivr, j
    2022
    2123      do i=ntlml,1,-1
    … …  
    4547      include 'INCL.H'
    4648
     49      integer i, ifirs, ilars, iat, j
     50
    4751      do i=1,ntlml
    4852

  • opereg.f

    r6650a56 r32289cd  
    44!
    55! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    6 !                      Shura Hayryan, Chin-Ku 
     6!                      Shura Hayryan, Chin-Ku
    77! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    88!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    3333      include 'INCL.H'
    3434      include 'INCP.H'
     35
     36      double precision xfvr, yfvr, zfvr, xfrvr, yfrvr, zfrvr, x1, y1, z1
     37      double precision a, sa, ca, xk, yk, zk, xi, yi, zi, xji, yji, zji
     38      double precision dx, dy, dz, x, y, z, xfat, yfat, zfat, xfrat
     39      double precision yfrat, zfrat, xb, yb, zb, ex, ey, ez, xfiv, yfiv
     40      double precision zfiv, xfriv, yfriv, zfriv
     41
     42      integer ntlvr, nml, ix2, ifivr, ilavr, i, ii, j, i1s, i1a, io, iv
     43      integer it, ia, ib, i2s, ims, i1, i2, i2a, iad, lad, ivw, i14
    3544
    3645      dimension xfat(mxat),yfat(mxat),zfat(mxat),
    … …  
    100109
    101110          eyrg = eyrg + xji**2 + yji**2 + zji**2  ! The regularization energy is just
    102                                                   ! the sum over the atom distances 
     111                                                  ! the sum over the atom distances
    103112                                                  ! squared.
    104113
    … …  
    303312      include 'INCL.H'
    304313
     314      double precision del, pro, gdn, enyreg
     315
     316      integer ii, nml, i
     317
    305318      parameter (del=1.d-7)
    306319
    … …  
    339352
    340353      include 'INCL.H'
     354
     355      double precision del, vlvrx, ovr, eynw, enyreg, gdn, gda
     356
     357      integer i, it, iv, nml
    341358
    342359      parameter (del=1.d-6)

  • opesol.f

    r6650a56 r32289cd  
    44!
    55! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    6 !                      Shura Hayryan, Chin-Ku 
     6!                      Shura Hayryan, Chin-Ku
    77! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    88!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    3939      integer ntlvr, nml, ix2, ifivr, ilavr, i, i1s, i1a, io, iv, it, ia
    4040      integer ib, i2s, ims, i1, i2, i2a, iad, lad, ivw, j, i14
    41      
     41
    4242      dimension xfat(mxat),yfat(mxat),zfat(mxat),
    4343     &          xfrat(mxat),yfrat(mxat),zfrat(mxat),
    … …  
    251251
    252252      include 'INCL.H'
    253      
     253
    254254      double precision del, vlvrx, ovr, eynw, esolan, gda, gdn
    255255

  • outpdb.f

    r6650a56 r32289cd  
    44!
    55! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    6 !                      Shura Hayryan, Chin-Ku 
     6!                      Shura Hayryan, Chin-Ku
    77! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    88!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    2424
    2525      include 'INCL.H'
     26
     27      double precision occ, bva
     28
     29      integer i0, i9, nml, im1, im2, ibegst, iout, iat, iml, irs, ifirs
     30      integer ifiat, nrs, i, j, iendst, nfi, ibd, jj, nbd
    2631
    2732      dimension ibd(4)
    … …  
    99104            call toupst(atnm)
    100105            call toupst(res)
    101  
     106
    102107            linout = ' '
    103108            write (linout,1) linty,iat,atnm,res(1:3),chid,irs,cdin,
    … …  
    170175      write (linout,'(a3)') 'END'
    171176      write(iout,'(a80)') linout
    172      
     177
    173178      close(iout)
    174179

  • outvar.f

    r6650a56 r32289cd  
    44!
    55! Copyright 2005       Frank Eisenmenger, U.H.E. Hansmann,
    6 !                      Shura Hayryan, Chin-Ku 
     6!                      Shura Hayryan, Chin-Ku
    77! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    88!                      Jan H. Meinke, Sandipan Mohanty
    99!
    1010! **************************************************************
    11  
     11
    1212      subroutine outvar(nml,fileName)
    13  
     13
    1414!--------------------------------------------------------------------
    15 !       Output of variables of the current protein conformation 
    16 !- 
     15!       Output of variables of the current protein conformation
     16!-
    1717!       nml != 0       :     molecule index
    1818!       nml  = 0       :     all molecules
    … …  
    2424
    2525      include 'INCL.H'
     26
     27      integer nml, ibegst, iout
    2628
    2729      character*(*) fileName
    … …  
    4951
    5052      include 'INCL.H'
     53
     54      integer nml, im1, im2, iml, i, iout, ibegst, ifivr, is, nursvr, it
    5155
    5256      character mlfd*7,strg(6)*17

  • partem_p.f

    r6650a56 r32289cd  
    4343
    4444      integer ifrrm, nmes, nswp, num_rep, i, j, nresi, iold, inode
    45       integer intem, iv, jold, idum1, idum2, idum3, mpi_integer
    46       integer mpi_comm_world, ierr, mpi_double_precision, nsw, nequi
     45      integer intem, iv, jold, idum1, idum2, idum3
     46      integer ierr, nsw, nequi
    4747      integer nml, nhel, mhel, nbet, mbet, mhb, imhb, nctot, ncnat
    48       integer mpi_comm_null, k1, k, nu, no1, in, jn
     48      integer k1, k, nu, no1, in, jn
    4949     
    5050      dimension  eavm(MAX_PROC),sph(MAX_PROC),intem(MAX_PROC),
    … …  
    225225            call flush(6)
    226226         endif
    227          if(mod(iold,nmes).eq.0) then
     227         if(mod(iold,nmes).eq.0.and.myrank.eq.0) then
    228228            if ((rep_id + 1).eq.trackID.and.myrank.eq.0) then
    229229               frame = iold /nmes

  • partem_s.f

    r6650a56 r32289cd  
    44!
    55! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    6 !                      Shura Hayryan, Chin-Ku 
     6!                      Shura Hayryan, Chin-Ku
    77! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    88!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    1010!
    1111! **************************************************************
    12      
     12
    1313      subroutine partem_s(num_rep, nequi, nswp, nmes, nsave, newsta,
    1414     &                     switch)
    … …  
    2020!
    2121      include 'INCL.H'
     22     
     23      double precision gasc, dv, grnd, vr, addang, energ, energy, acz
     24      double precision ee, rgy, temp0, delta, ra
     25
     26      integer i, j, nml, k, nstart, iv, nsw, iswitch, k1, num_rep1, nu
     27      integer in, jn
     28
    2229      external can_weight
    2330! TODO Store global coordinates in pgbpr
    … …  
    4451!     ipoi:   Points to replica ipoi(k) which is currently
    4552!             at inverse temperature pbe(k)
    46 !     eol:    energy of each replica 
     53!     eol:    energy of each replica
    4754!     acc:    accepatance rate of each replica
    4855!
    … …  
    5562      allocate(temp(num_rep), pbe(num_rep), eol(num_rep))
    5663      allocate(acc(num_rep), ipoi(num_rep))
    57      
     64
    5865      if (.not.(allocated(coor_G).and.allocated(temp)
    5966     &          .and. allocated(pbe).and. allocated(eol)
    … …  
    7784         ipoi(i) = 0
    7885      end do
    79        
     86
    8087
    8188!     READ IN TEMPERATURES
    … …  
    9097
    9198      open(13,file='time.d',status='unknown')
    92      
     99
    93100      if(.not.newsta) then
    94101! READ START Values
    … …  
    122129                  pgbpr(j, nml, i) = gbpr(j, nml)
    123130               end do
    124             end do           
    125          end do
    126      
     131            end do
     132         end do
     133
    127134!      Equilibrization for each replica (No replica exchange move)
    128135         do k=1,num_rep
    … …  
    141148               CALL METROPOLIS(energ,acz,can_weight)
    142149            end do
    143             write(*,*) 'Start energy after equilibration for replica:', 
     150            write(*,*) 'Start energy after equilibration for replica:',
    144151     &                 k, energ
    145152            do i=1,nvr
    … …  
    155162         end do
    156163      end if
    157            
     164
    158165! Now begins the simulation with Multiple Markov Chains
    159166      iswitch = 1
    … …  
    176183            end do
    177184            CALL METROPOLIS(energ,acz,can_weight)
    178 ! 
     185!
    179186            if(mod(nsw,nmes).eq.0) then
    180187! Measure and store here all quantities you want to analyse later
    … …  
    195202               end do
    196203            end do
    197        
     204
    198205            eol(k) = energ
    199206            acc(k) = acz
    … …  
    214221            j=i+1
    215222            if(i.eq.num_rep) j=1
    216             in=ipoi(i)     
     223            in=ipoi(i)
    217224            jn=ipoi(j)
    218225            delta=-pbe(i)*eol(jn)-pbe(j)*eol(in)

  • pdbread.f

    r6650a56 r32289cd  
    44!
    55! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    6 !                      Shura Hayryan, Chin-Ku 
     6!                      Shura Hayryan, Chin-Ku
    77! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    88!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    2121! ......................................................
    2222
    23       implicit real*8 (a-h,o-z)
    24       implicit integer*4 (i-n)
    25 
    2623      include 'INCP.H'
     24     
     25      double precision cor
     26
     27      integer ier, l, iendst, lunpdb, io, iopfil, iat, i
    2728
    2829! -------------------------- input
    … …  
    4748        lunpdb = 99
    4849      else
    49         write (*,'(a)') 
     50        write (*,'(a)')
    5051     &    ' pdbread> empty file name to read pdb-structure'
    5152
    … …  
    5657
    5758      if (io.le.0) then
    58         write (*,'(a,/,a)') 
     59        write (*,'(a,/,a)')
    5960     &  ' pdbread> ERROR opening file to read pdb-structure: ',
    6061     &  pdbfil(1:iendst(pdbfil))
    … …  
    6970
    7071      if ( line(17:17).ne.' ' )  then
    71         write (*,'(a,/,a,/,a,/,2a)') 
     72        write (*,'(a,/,a,/,a,/,2a)')
    7273     &  ' pdbread> found alternate atom location: ',
    7374     &  '                !',
    … …  
    8586
    8687      if ((natp+1).gt.MXATP) then
    87         write (*,'(a,i5,a,/,a)') 
     88        write (*,'(a,i5,a,/,a)')
    8889     &  ' pdbread>  >MXATP (',MXATP,') ATOM lines in PDB file ',
    8990     &  pdbfil(1:iendst(pdbfil))
    … …  
    9697
    9798        if ((nchp+1).gt.MXCHP) then
    98           write (*,'(a,i3,a,/,a)') 
     99          write (*,'(a,i3,a,/,a)')
    99100     &    ' pdbread>  >MXCHP (',MXCHP,') chains in PDB file ',
    100101     &    pdbfil(1:iendst(pdbfil))
    … …  
    105106
    106107        if ((nrsp+1).gt.MXRSP) then
    107           write (*,'(a,i3,a,/,a)') 
     108          write (*,'(a,i3,a,/,a)')
    108109     &    ' pdbread>  >MXRSP (',MXRSP,') residues in PDB file ',
    109110     &    pdbfil(1:iendst(pdbfil))
    … …  
    140141
    141142        if ((nrsp+1).gt.MXRSP) then
    142           write (*,'(a,i3,a,/,a)') 
     143          write (*,'(a,i3,a,/,a)')
    143144     &    ' pdbread>  >MXRSP (',MXRSP,') residues in PDB file ',
    144145     &    pdbfil(1:iendst(pdbfil))
    … …  
    171172      goto 1
    172173
    173     2 write (*,'(a,/,a,/,2a)') 
     174    2 write (*,'(a,/,a,/,2a)')
    174175     &  ' pdbread> ERROR reading ATOM line ',
    175176     &  line(:l),
    … …  
    194195      else
    195196
    196         write (*,'(a,/,a)') 
     197        write (*,'(a,/,a)')
    197198     &  ' pdbread> NO atom coordinates selected from file ',
    198199     &  pdbfil(1:iendst(pdbfil))
    … …  
    226227      include 'INCL.H'
    227228      include 'INCP.H'
     229
     230      double precision dihedr, rm, av1, av2, rmsd, h, x, y, z
     231
     232      integer nml, nrs, nc, irb, ire, irs, i, ii, it, i1, i2, i3, i4, j1
     233      integer j2, j3, j4, inew, j, ix
    228234
    229235      character res*3
    … …  
    451457      include 'INCP.H'
    452458
     459      integer irs, nml, i1, i2, iat, ix, i
     460
    453461      character*4 atm
    454462
    … …  
    476484              endif
    477485            enddo
    478                
     486
    479487!            write(*,'(8a)') ' pdbvars> ',atm,' not found in '
    480488!     #       ,chnp(nc),' ',rsidp(irs),' ',rsnmp(irs)
    … …  
    484492    1     ixatp(iat)=ix
    485493
    486         enddo   ! SMMP atoms of 'irs' 
     494        enddo   ! SMMP atoms of 'irs'
    487495      enddo   ! residues
    488496
    … …  
    494502      include 'INCL.H'
    495503
     504      double precision tol, h1, h2, h3, x1, x2, x3, d, z1, z2, z3, y1
     505      double precision y2, y3, yp1, yp2
     506
     507      integer i1, nml, i2, i3, i
     508
    496509      parameter (TOL = 1.d-12)
    497510
    … …  
    499512! -> determine global parameters: shifts dX,dY,dZ
    500513! & angles alpha,beta,gamma [rad], put into 'gbpr'
    501 ! 
     514!
    502515! CALLS: none
    503516!
    … …  
    547560      y2=z3*x1-z1*x3
    548561      y3=z1*x2-z2*x1
    549      
     562
    550563      if ( ( 1.d0 - abs(y3) ) .gt. TOL )  then       ! ============ |beta| < PI/2
    551564

  • regul.f

    r6650a56 r32289cd  
    44!
    55! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    6 !                      Shura Hayryan, Chin-Ku 
     6!                      Shura Hayryan, Chin-Ku
    77! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    88!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    1010! **************************************************************
    1111
    12      
     12
    1313      subroutine regul(nml, iter, nsteps, acc)
    1414
    … …  
    2626      include 'INCL.H'
    2727      include 'INCP.H'
    28      
     28
    2929!f2py intent(in) nml
    3030!f2py intent(in) iter
    3131!f2py intent(in) nsteps
    3232!f2py intent(in) acc
     33
     34      double precision acc, rm, av1, av2, rmsd, dn
     35
     36      integer nsteps, nml, nrs, i, n, iter, it
    3337
    3438      dimension rm(3,3),av1(3),av2(3)
    … …  
    6064
    6165
    62       do i = ivrml1(nml),nvrml(nml) 
     66      do i = ivrml1(nml),nvrml(nml)
    6367        fxvro(i) = fxvr(i)  ! save
    6468        if (isrfvr(i)) fxvr(i) = .true.  ! fix vars. defined in ref.str.
    … …  
    7781      call cnteny(nml)
    7882
    79       do i = ivrml1(nml),nvrml(nml) 
     83      do i = ivrml1(nml),nvrml(nml)
    8084        fxvr(i) = fxvro(i)  ! restore
    8185      enddo  ! vars.
    … …  
    8690      wtey = 0.d0
    8791
    88       n=iter           
     92      n=iter
    8993      dn=1.d0/dble(n)
    9094

  • rgyr.f

    r6650a56 r32289cd  
    2525!
    2626      include 'INCL.H'
     27     
     28      double precision dn, dnp, dnh, dx, dxp, dy, dyp, dyh, dz, dzp, dzh
     29      double precision d2, d2p, d2h, xi, yi, zi, dxh, rg2, rg2p, rg2h
     30      double precision rgy, ee
     31
     32      integer nml, nml1, nml2, nat, i, i1, i2
     33
    2734!f2py intent(in) nml
    2835!f2py intent(out) rgy

  • rmsdfun.f

    r6650a56 r32289cd  
    55!
    66! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
    7 !                      Shura Hayryan, Chin-Ku 
     7!                      Shura Hayryan, Chin-Ku
    88! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
    99!                      Jan H. Meinke, Sandipan Mohanty
    … …  
    2626!
    2727! Input
     28      double precision xrf, yrf, zrf, rm, av1, av2, rssd
     29
     30      integer nml, ir1, ir2, ixat, isl
     31
    2832      dimension ixat(mxat),xrf(mxatp),yrf(mxatp),zrf(mxatp)
    2933! Local
    … …  
    5761!
    5862!      NB  uncomment last lines in 'fitmol' to return coordinates
    59 !          in 'x2' after fitting the ref. str. onto SMMP structure 
     63!          in 'x2' after fitting the ref. str. onto SMMP structure
    6064! ----------------------------------------------------------------
    6165
    6266      include 'INCL.H'
    6367      include 'INCP.H'
     68     
     69      double precision x1, x2, xrf, yrf, zrf, rm, av1, av2, rmsd
     70
     71      integer nml, nr, na, n, im, ir, ia, ir1, ir2, ix, ixat, isl
     72
    6473
    6574!-------------------------------------------------------- input
    … …  
    147156! .......................................................
    148157!f2py intent(out) rmsd
    149                    
     158
    150159      include 'INCL.H'
    151160!      implicit real*8 (a-h,o-z)
    152161!      implicit integer*4 (i-n)
    153  
     162
    154163! ------------------------------------------- input/output
     164      double precision dn, a1, a2, x1, x2, q, dm, dp, dxm, dym, dzm, dxp
     165      double precision dyp, dzp, e, v, em, rmsd, rm
     166
     167      integer n, i, j, ndim4, im
     168
    155169      dimension x1(3,mxat),x2(3,mxat)
    156170! -------------------------------------------------- local
    … …  
    270284!f2py intent(out) d
    271285!f2py intent(out) v
     286      integer nmax
     287
    272288      parameter (NMAX=500)
    273      
     289
    274290      integer n,nrot,i,ip,iq,j
    275291
    … …  
    398414!
    399415!------------------------------------------------------------------------------
    400 ! Reads in pdb-file 'string' into INCP.H and initalizes 
     416! Reads in pdb-file 'string' into INCP.H and initalizes
    401417! the files that 'rmdsopt' needs to calculate the rmsd
    402418! of a configuration with the pdb-configuration
    … …  
    409425      include 'INCP.H'
    410426
     427      integer i, nml, ier
     428
    411429      character string*(*)
    412  
     430
    413431      if(string.eq.'smmp') then
    414432!
    … …  
    428446!
    429447         call pdbread(string,ier)
    430          if(ier.ne.0) stop 
     448         if(ier.ne.0) stop
    431449         call atixpdb(nml)
    432450!

  • twister.f

    r6650a56 r32289cd  
    4747* May 03, 2000
    4848*
    49 * compile with: f90 -i32 -dp -o mt19937 mt19937.f 
     49* compile with: f90 -i32 -dp -o mt19937 mt19937.f
    5050* Changes:
    5151*    Do not use UMASK as parameter due to restrictions of CF90
    … …  
    5555*
    5656*    Now the T3E produces the same random number as a RS6000 running AIX
    57 *     
     57*
    5858* from other programs sgrnd() must be called with integer*4
    5959*
    … …  
    6464      subroutine sgrnd(seed)
    6565*
    66       implicit integer(a-z)
     66      integer seed, tmaskb, tmaskc, tshftu, tshfts, tshftt
     67      integer tshftl, y
     68
     69      integer N, N1, mti, mt, m, mata, lmask, mag01
     70      integer kk
    6771*
    6872* Period parameters
    … …  
    8892
    8993      block data twbloks
     94      integer n,n1,mt,mti
    9095      parameter(N     =  624)
    9196      parameter(N1    =  N+1)
    … …  
    99104      double precision function grnd()
    100105*
    101       implicit integer(a-z)
     106      integer tmaskb, tmaskc, tshftu, tshfts, tshftt, tshftl, y
     107
     108      integer n, n1, mti, mt, m, mata, lmask, mag01, kk
     109
    102110*
    103111* Period parameters

  • utilities.f

    r6650a56 r32289cd  
    88
    99!! Calculate the best way to distribute the work load across processors.
    10 !  It calculates the average number of interactions and then tries to 
     10!  It calculates the average number of interactions and then tries to
    1111!  assign a number of interactions to each processor that is as close
    1212!  as possible to the average. The routine should be called once for
    … …  
    2121      include 'INCL.H'
    2222
     23      integer i1ms, io, iv, i2ms, ms
     24
    2325      integer num_ppr, nml
    2426      integer idxOfFirstVariable, idxOfLastVariable
    … …  
    2628      integer totalct, irank, itarget
    2729      double precision ipps
    28      
     30
    2931      if (nml.eq.0) then
    3032         idxOfFirstVariable = ivrml1(1)
    … …  
    3638            end do
    3739         end do
    38       else 
     40      else
    3941         idxOfFirstVariable = ivrml1(nml)
    4042         idxOfLastVariable = ivrml1(nml) + nvrml(nml) - 1
    … …  
    4446         end do
    4547      end if
    46      
     48
    4749      isum = 0
    4850      do io = idxOfLastVariable, idxOfFirstVariable, - 1
    … …  
    5254         do ms = i1ms, i2ms
    5355            do at = latms1(ms), latms2(ms)
    54                do ivw=ivwat1(at),ivwat2(at) 
     56               do ivw=ivwat1(at),ivwat2(at)
    5557                  do j=lvwat1(ivw),lvwat2(ivw)
    5658                     isum = isum + 1
    5759                  end do
    5860               end do
    59                do i14=i14at1(at),i14at2(at)   
     61               do i14=i14at1(at),i14at2(at)
    6062                  isum = isum + 1
    6163               end do
    … …  
    6668      write (*,*) "Total number of interactions:", isum
    6769      write (*,*) "Average # of interactions per processor", ipps
    68      
     70
    6971      totalct = 0
    7072      irank = 1
    … …  
    7274      if (nml.eq.0) then
    7375         i1ms = imsml1(ntlml)+ nmsml(ntlml)
    74       else 
     76      else
    7577         i1ms = imsml1(nml)+ nmsml(nml)
    7678      end if
    7779      do io = idxOfLastVariable, idxOfFirstVariable, - 1
    78          isum = 0 
     80         isum = 0
    7981         iv = iorvr(io)
    8082         i2ms = i1ms - 1
    … …  
    8385            do at = latms1(ms), latms2(ms)
    8486               atct = atct + 1
    85                do ivw=ivwat1(at),ivwat2(at) 
     87               do ivw=ivwat1(at),ivwat2(at)
    8688                  do j=lvwat1(ivw),lvwat2(ivw)
    8789                     isum = isum + 1
    8890                  end do
    8991               end do
    90                do i14=i14at1(at),i14at2(at)   
     92               do i14=i14at1(at),i14at2(at)
    9193                  isum = isum + 1
    9294               end do
    … …  
    113115
    114116      end subroutine distributeWorkLoad
    115      
     117
    116118!-----------------------------------------------------------------------
    117119!     The function fileNameMP takes a template of a file name in the
    … …  
    128130!     \endcode
    129131!     will output base_0011.dat.
    130 !     
     132!
    131133!     @param base the base file name, e.g., base_0000.dat.
    132134!     @param i1 index of the first character that may be replaced
    133135!     @param i2 index of the last character that may be replaced
    134136!     @param rank the number that should be inserted into the file name.
    135 !     
     137!
    136138!     @return file name for rank
    137139!-----------------------------------------------------------------------
    … …  
    165167         write(fileNameMP(i2-5:i2), '(i6)') rank
    166168      endif
    167       end function fileNameMP 
     169      end function fileNameMP
    168170!     End fileNameMP
    169171
    … …  
    172174!     Add messages to log. This routine takes the log (debugging) mes-
    173175!     sages and writes them to the log file if the log level is less or
    174 !     equal to the maximum log level given by the global variable 
     176!     equal to the maximum log level given by the global variable
    175177!     MAXLOGLEVEL.
    176178!
    177179!     @author Jan H. Meinke
    178180!
    179 !     @param loglevel level at which this message should be added to 
     181!     @param loglevel level at which this message should be added to
    180182!            the log.
    181183!     @param message message to be written to the log.
    182 !     @param rank global rank of this node if running an MPI job zero 
     184!     @param rank global rank of this node if running an MPI job zero
    183185!            otherwise.
    184186!----------------------------------------------------------------------
    185187      subroutine addLogMessage(loglevel, message, rank)
    186      
     188
     189      integer maxloglevel, logfileunit
     190
    187191         integer :: loglevel, rank
    188192         character(LEN=*) :: message
    189          
     193
    190194         if (loglevel <= MAXLOGLEVEL) then
    191195            write(LOGFILEUNIT, *) message
    192196         end if
    193          
     197
    194198      end subroutine addLogMessage

  • zimmer.f

    r6650a56 r32289cd  
    4040!f2py intent(in) nresi
    4141      character*1 zim
     42      double precision xphi, xpsi
     43
     44      integer j, nresi, i, iv, nres, nursvr
    4245
    4346     
Note: See TracChangeset for help on using the changeset viewer.