source/wham/src-restraints/parmread.F

   1       subroutine parmread(iparm,*)
   2 C
   3 C Read the parameters of the probability distributions of the virtual-bond
   4 C valence angles and the side chains and energy parameters.
   5 C
   6       implicit real*8 (a-h,o-z)
   7       include 'DIMENSIONS'
   8       include 'DIMENSIONS.ZSCOPT'
   9       include 'DIMENSIONS.FREE'
  10       include 'COMMON.IOUNITS'
  11       include 'COMMON.CHAIN'
  12       include 'COMMON.INTERACT'
  13       include 'COMMON.GEO'
  14       include 'COMMON.LOCAL'
  15       include 'COMMON.TORSION'
  16       include 'COMMON.FFIELD'
  17       include 'COMMON.NAMES'
  18       include 'COMMON.SBRIDGE'
  19       include 'COMMON.WEIGHTS'
  20       include 'COMMON.ENEPS'
  21       include 'COMMON.SCCOR'
  22       include 'COMMON.SCROT'
  23       include 'COMMON.FREE'
  24       character*1 t1,t2,t3
  25       character*1 onelett(4) /"G","A","P","D"/
  26       logical lprint
  27       dimension blower(3,3,maxlob)
  28       character*800 controlcard
  29       character*256 bondname_t,thetname_t,rotname_t,torname_t,
  30      & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
  31      & sccorname_t
  32       integer ilen
  33       external ilen
  34       character*16 key
  35       integer iparm
  36       double precision ip,mp
  37 C
  38 C Body
  39 C
  40 C Set LPRINT=.TRUE. for debugging
  41       dwa16=2.0d0**(1.0d0/6.0d0)
  42       lprint=.false.
  43       itypro=20
  44 C Assign virtual-bond length
  45       vbl=3.8D0
  46       vblinv=1.0D0/vbl
  47       vblinv2=vblinv*vblinv
  48       call card_concat(controlcard,.true.)
  49       wname(4)="WCORRH"
  50       do i=1,n_ene
  51         key = wname(i)(:ilen(wname(i)))
  52         call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
  53       enddo
  54       call reada(controlcard,"D0CM",d0cm,3.78d0)
  55       call reada(controlcard,"AKCM",akcm,15.1d0)
  56       call reada(controlcard,"AKTH",akth,11.0d0)
  57       call reada(controlcard,"AKCT",akct,12.0d0)
  58       call reada(controlcard,"V1SS",v1ss,-1.08d0)
  59       call reada(controlcard,"V2SS",v2ss,7.61d0)
  60       call reada(controlcard,"V3SS",v3ss,13.7d0)
  61       call reada(controlcard,"EBR",ebr,-5.50D0)
  62 c      dyn_ss=(index(controlcard,'DYN_SS').gt.0)
  63       write (iout,*) "iparm",iparm," myparm",myparm
  64 c      do i=1,maxres
  65 c        dyn_ss_mask(i)=.false.
  66 c      enddo
  67       do i=1,maxres-1
  68         do j=i+1,maxres
  69           dyn_ssbond_ij(i,j)=1.0d300
  70         enddo
  71       enddo
  72       call reada(controlcard,"HT",Ht,0.0D0)
  73
  74 c      if(me.eq.king.or..not.out1file) then
  75 c       print *,'indpdb=',indpdb,' pdbref=',pdbref
  76 c      endif
  77 c If reading not own parameters, skip assignment
  78 cc      write(iout,*) "KURWA", ww(15)
  79       if (iparm.eq.myparm .or. .not.separate_parset) then
  80
  81 c
  82 c Setup weights for UNRES
  83 c
  84       wsc=ww(1)
  85       wscp=ww(2)
  86       welec=ww(3)
  87       wcorr=ww(4)
  88       wcorr5=ww(5)
  89       wcorr6=ww(6)
  90       wel_loc=ww(7)
  91       wturn3=ww(8)
  92       wturn4=ww(9)
  93       wturn6=ww(10)
  94       wang=ww(11)
  95       wscloc=ww(12)
  96       wtor=ww(13)
  97       wtor_d=ww(14)
  98       wvdwpp=ww(16)
  99       wstrain=ww(15)
 100       wbond=ww(18)
 101       wsccor=ww(19)
 102
 103       endif
 104 cc      write(iout,*) "KURWA", wstrain,akcm,akth,wsc,dyn_ss
 105
 106       call card_concat(controlcard,.false.)
 107
 108 c Return if not own parameters
 109
 110       if (iparm.ne.myparm .and. separate_parset) return
 111
 112       call reads(controlcard,"BONDPAR",bondname_t,bondname)
 113       open (ibond,file=bondname_t,status='old')
 114       rewind(ibond)
 115       call reads(controlcard,"THETPAR",thetname_t,thetname)
 116       open (ithep,file=thetname_t,status='old')
 117       rewind(ithep)
 118       call reads(controlcard,"ROTPAR",rotname_t,rotname)
 119       open (irotam,file=rotname_t,status='old')
 120       rewind(irotam)
 121       call reads(controlcard,"TORPAR",torname_t,torname)
 122       open (itorp,file=torname_t,status='old')
 123       rewind(itorp)
 124       call reads(controlcard,"TORDPAR",tordname_t,tordname)
 125       open (itordp,file=tordname_t,status='old')
 126       rewind(itordp)
 127       call reads(controlcard,"SCCORAR",sccorname_t,sccorname)
 128       open (isccor,file=sccorname_t,status='old')
 129       rewind(isccor)
 130       call reads(controlcard,"FOURIER",fouriername_t,fouriername)
 131       open (ifourier,file=fouriername_t,status='old')
 132       rewind(ifourier)
 133       call reads(controlcard,"ELEPAR",elename_t,elename)
 134       open (ielep,file=elename_t,status='old')
 135       rewind(ielep)
 136       call reads(controlcard,"SIDEPAR",sidename_t,sidename)
 137       open (isidep,file=sidename_t,status='old')
 138       rewind(isidep)
 139       call reads(controlcard,"SCPPAR",scpname_t,scpname)
 140       open (iscpp,file=scpname_t,status='old')
 141       rewind(iscpp)
 142       write (iout,*) "Parameter set:",iparm
 143       write (iout,*) "Energy-term weights:"
 144       do i=1,n_ene
 145         write (iout,'(a16,f10.5)') wname(i),ww(i)
 146       enddo
 147       write (iout,*) "Sidechain potential file        : ",
 148      &  sidename_t(:ilen(sidename_t))
 149 #ifndef OLDSCP
 150       write (iout,*) "SCp potential file              : ",
 151      &  scpname_t(:ilen(scpname_t))
 152 #endif
 153       write (iout,*) "Electrostatic potential file    : ",
 154      &  elename_t(:ilen(elename_t))
 155       write (iout,*) "Cumulant coefficient file       : ",
 156      &  fouriername_t(:ilen(fouriername_t))
 157       write (iout,*) "Torsional parameter file        : ",
 158      &  torname_t(:ilen(torname_t))
 159       write (iout,*) "Double torsional parameter file : ",
 160      &  tordname_t(:ilen(tordname_t))
 161       write (iout,*) "Backbone-rotamer parameter file : ",
 162      &  sccorname(:ilen(sccorname))
 163       write (iout,*) "Bond & inertia constant file    : ",
 164      &  bondname_t(:ilen(bondname_t))
 165       write (iout,*) "Bending parameter file          : ",
 166      &  thetname_t(:ilen(thetname_t))
 167       write (iout,*) "Rotamer parameter file          : ",
 168      &  rotname_t(:ilen(rotname_t))
 169
 170 c
 171 c Read the virtual-bond parameters, masses, and moments of inertia
 172 c and Stokes' radii of the peptide group and side chains
 173 c
 174 #ifdef CRYST_BOND
 175       read (ibond,*) vbldp0,akp
 176       do i=1,ntyp
 177         nbondterm(i)=1
 178         read (ibond,*) vbldsc0(1,i),aksc(1,i)
 179         dsc(i) = vbldsc0(1,i)
 180         if (i.eq.10) then
 181           dsc_inv(i)=0.0D0
 182         else
 183           dsc_inv(i)=1.0D0/dsc(i)
 184         endif
 185       enddo
 186 #else
 187       read (ibond,*) ijunk,vbldp0,akp,rjunk
 188       do i=1,ntyp
 189         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
 190      &   j=1,nbondterm(i))
 191         dsc(i) = vbldsc0(1,i)
 192         if (i.eq.10) then
 193           dsc_inv(i)=0.0D0
 194         else
 195           dsc_inv(i)=1.0D0/dsc(i)
 196         endif
 197       enddo
 198 #endif
 199       if (lprint) then
 200         write(iout,'(/a/)')"Force constants virtual bonds:"
 201         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
 202      &   'inertia','Pstok'
 203         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
 204         do i=1,ntyp
 205           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
 206      &      vbldsc0(1,i),aksc(1,i),abond0(1,i)
 207           do j=2,nbondterm(i)
 208             write (iout,'(13x,3f10.5)')
 209      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
 210           enddo
 211         enddo
 212       endif
 213 #ifdef CRYST_THETA
 214 C
 215 C Read the parameters of the probability distribution/energy expression
 216 C of the virtual-bond valence angles theta
 217 C
 218       do i=1,ntyp
 219         read (ithep,*) a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
 220         read (ithep,*) (polthet(j,i),j=0,3)
 221         read (ithep,*) (gthet(j,i),j=1,3)
 222         read (ithep,*) theta0(i),sig0(i),sigc0(i)
 223         sigc0(i)=sigc0(i)**2
 224       enddo
 225       close (ithep)
 226       if (lprint) then
 227 c       write (iout,'(a)')
 228 c    &    'Parameters of the virtual-bond valence angles:'
 229 c       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
 230 c    & '    ATHETA0   ','         A1   ','        A2    ',
 231 c    & '        B1    ','         B2   '
 232 c       do i=1,ntyp
 233 c         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 234 c    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
 235 c       enddo
 236 c       write (iout,'(/a/9x,5a/79(1h-))')
 237 c    & 'Parameters of the expression for sigma(theta_c):',
 238 c    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
 239 c    & '     ALPH3    ','    SIGMA0C   '
 240 c       do i=1,ntyp
 241 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 242 c    &      (polthet(j,i),j=0,3),sigc0(i)
 243 c       enddo
 244 c       write (iout,'(/a/9x,5a/79(1h-))')
 245 c    & 'Parameters of the second gaussian:',
 246 c    & '    THETA0    ','     SIGMA0   ','        G1    ',
 247 c    & '        G2    ','         G3   '
 248 c       do i=1,ntyp
 249 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
 250 c    &       sig0(i),(gthet(j,i),j=1,3)
 251 c       enddo
 252         write (iout,'(a)')
 253      &    'Parameters of the virtual-bond valence angles:'
 254         write (iout,'(/a/9x,5a/79(1h-))')
 255      & 'Coefficients of expansion',
 256      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
 257      & '   b1*10^1    ','    b2*10^1   '
 258         do i=1,ntyp
 259           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
 260      &        a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
 261         enddo
 262         write (iout,'(/a/9x,5a/79(1h-))')
 263      & 'Parameters of the expression for sigma(theta_c):',
 264      & ' alpha0       ','  alph1       ',' alph2        ',
 265      & ' alhp3        ','   sigma0c    '
 266         do i=1,ntyp
 267           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
 268      &      (polthet(j,i),j=0,3),sigc0(i)
 269         enddo
 270         write (iout,'(/a/9x,5a/79(1h-))')
 271      & 'Parameters of the second gaussian:',
 272      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
 273      & '        G2    ','   G3*10^1    '
 274         do i=1,ntyp
 275           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
 276      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 277         enddo
 278       endif
 279 #else
 280 C
 281 C Read the parameters of Utheta determined from ab initio surfaces
 282 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 283 C
 284       read (ithep,*) nthetyp,ntheterm,ntheterm2,
 285      &  ntheterm3,nsingle,ndouble
 286       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 287       read (ithep,*) (ithetyp(i),i=1,ntyp1)
 288       do i=1,maxthetyp
 289         do j=1,maxthetyp
 290           do k=1,maxthetyp
 291             aa0thet(i,j,k)=0.0d0
 292             do l=1,ntheterm
 293               aathet(l,i,j,k)=0.0d0
 294             enddo
 295             do l=1,ntheterm2
 296               do m=1,nsingle
 297                 bbthet(m,l,i,j,k)=0.0d0
 298                 ccthet(m,l,i,j,k)=0.0d0
 299                 ddthet(m,l,i,j,k)=0.0d0
 300                 eethet(m,l,i,j,k)=0.0d0
 301               enddo
 302             enddo
 303             do l=1,ntheterm3
 304               do m=1,ndouble
 305                 do mm=1,ndouble
 306                  ffthet(mm,m,l,i,j,k)=0.0d0
 307                  ggthet(mm,m,l,i,j,k)=0.0d0
 308                 enddo
 309               enddo
 310             enddo
 311           enddo
 312         enddo
 313       enddo
 314       do i=1,nthetyp
 315         do j=1,nthetyp
 316           do k=1,nthetyp
 317             read (ithep,'(3a)') res1,res2,res3
 318             read (ithep,*) aa0thet(i,j,k)
 319             read (ithep,*)(aathet(l,i,j,k),l=1,ntheterm)
 320             read (ithep,*)
 321      &       ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
 322      &        (ccthet(lll,ll,i,j,k),lll=1,nsingle),
 323      &        (ddthet(lll,ll,i,j,k),lll=1,nsingle),
 324      &        (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
 325             read (ithep,*)
 326      &      (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
 327      &         ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
 328      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 329           enddo
 330         enddo
 331       enddo
 332 C
 333 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 334 C coefficients of theta-and-gamma-dependent terms are zero.
 335 C
 336       do i=1,nthetyp
 337         do j=1,nthetyp
 338           do l=1,ntheterm
 339             aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
 340             aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
 341           enddo
 342           aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
 343           aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
 344         enddo
 345         do l=1,ntheterm
 346           aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
 347         enddo
 348         aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
 349       enddo
 350 C
 351 C Control printout of the coefficients of virtual-bond-angle potentials
 352 C
 353       if (lprint) then
 354         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 355         do i=1,nthetyp+1
 356           do j=1,nthetyp+1
 357             do k=1,nthetyp+1
 358               write (iout,'(//4a)')
 359      &         'Type ',onelett(i),onelett(j),onelett(k)
 360               write (iout,'(//a,10x,a)') " l","a[l]"
 361               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
 362               write (iout,'(i2,1pe15.5)')
 363      &           (l,aathet(l,i,j,k),l=1,ntheterm)
 364             do l=1,ntheterm2
 365               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
 366      &          "b",l,"c",l,"d",l,"e",l
 367               do m=1,nsingle
 368                 write (iout,'(i2,4(1pe15.5))') m,
 369      &          bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
 370      &          ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
 371               enddo
 372             enddo
 373             do l=1,ntheterm3
 374               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
 375      &          "f+",l,"f-",l,"g+",l,"g-",l
 376               do m=2,ndouble
 377                 do n=1,m-1
 378                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
 379      &              ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
 380      &              ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
 381                 enddo
 382               enddo
 383             enddo
 384           enddo
 385         enddo
 386       enddo
 387       call flush(iout)
 388       endif
 389 #endif
 390
 391 #ifdef CRYST_SC
 392 C
 393 C Read the parameters of the probability distribution/energy expression
 394 C of the side chains.
 395 C
 396       do i=1,ntyp
 397         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
 398         if (i.eq.10) then
 399           dsc_inv(i)=0.0D0
 400         else
 401           dsc_inv(i)=1.0D0/dsc(i)
 402         endif
 403         if (i.ne.10) then
 404         do j=1,nlob(i)
 405           do k=1,3
 406             do l=1,3
 407               blower(l,k,j)=0.0D0
 408             enddo
 409           enddo
 410         enddo
 411         bsc(1,i)=0.0D0
 412         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
 413         do j=2,nlob(i)
 414           read (irotam,*) bsc(j,i)
 415           read (irotam,*) (censc(k,j,i),k=1,3),
 416      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 417         enddo
 418         do j=1,nlob(i)
 419           do k=1,3
 420             do l=1,k
 421               akl=0.0D0
 422               do m=1,3
 423                 akl=akl+blower(k,m,j)*blower(l,m,j)
 424               enddo
 425               gaussc(k,l,j,i)=akl
 426               gaussc(l,k,j,i)=akl
 427             enddo
 428           enddo
 429         enddo
 430         endif
 431       enddo
 432       close (irotam)
 433       if (lprint) then
 434         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 435         do i=1,ntyp
 436           nlobi=nlob(i)
 437           if (nlobi.gt.0) then
 438           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 439      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 440 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 441 c          write (iout,'(a,f10.4,4(16x,f10.4))')
 442 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
 443 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
 444            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 445      &                             'log h',(bsc(j,i),j=1,nlobi)
 446            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 447      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 448 c          write (iout,'(a)')
 449 c         do j=1,nlobi
 450 c           ind=0
 451 c           do k=1,3
 452 c             do l=1,k
 453 c              ind=ind+1
 454 c              blower(k,l,j)=gaussc(ind,j,i)
 455 c             enddo
 456 c           enddo
 457 c         enddo
 458           do k=1,3
 459             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 460      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 461           enddo
 462           endif
 463         enddo
 464       endif
 465 #else
 466 C
 467 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 468 C added by Urszula Kozlowska 07/11/2007
 469 C
 470       do i=1,ntyp
 471         read (irotam,*)
 472        if (i.eq.10) then
 473          read (irotam,*)
 474        else
 475          do j=1,65
 476            read(irotam,*) sc_parmin(j,i)
 477          enddo
 478        endif
 479       enddo
 480 #endif
 481       close(irotam)
 482 #ifdef CRYST_TOR
 483 C
 484 C Read torsional parameters in old format
 485 C
 486       read (itorp,*) ntortyp,nterm_old
 487       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 488       read (itorp,*) (itortyp(i),i=1,ntyp)
 489       do i=1,ntortyp
 490         do j=1,ntortyp
 491           read (itorp,'(a)')
 492           do k=1,nterm_old
 493             read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
 494           enddo
 495         enddo
 496       enddo
 497       close (itorp)
 498       if (lprint) then
 499         write (iout,'(/a/)') 'Torsional constants:'
 500         do i=1,ntortyp
 501           do j=1,ntortyp
 502             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
 503             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
 504           enddo
 505         enddo
 506       endif
 507
 508
 509 #else
 510 C
 511 C Read torsional parameters
 512 C
 513       read (itorp,*) ntortyp
 514       read (itorp,*) (itortyp(i),i=1,ntyp)
 515       write (iout,*) 'ntortyp',ntortyp
 516       do i=1,ntortyp
 517         do j=1,ntortyp
 518           read (itorp,*) nterm(i,j),nlor(i,j)
 519           v0ij=0.0d0
 520           si=-1.0d0
 521           do k=1,nterm(i,j)
 522             read (itorp,*) kk,v1(k,i,j),v2(k,i,j)
 523             v0ij=v0ij+si*v1(k,i,j)
 524             si=-si
 525           enddo
 526           do k=1,nlor(i,j)
 527             read (itorp,*) kk,vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 528             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
 529           enddo
 530           v0(i,j)=v0ij
 531         enddo
 532       enddo
 533       close (itorp)
 534       if (lprint) then
 535         write (iout,'(/a/)') 'Torsional constants:'
 536         do i=1,ntortyp
 537           do j=1,ntortyp
 538             write (iout,*) 'ityp',i,' jtyp',j
 539             write (iout,*) 'Fourier constants'
 540             do k=1,nterm(i,j)
 541               write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
 542             enddo
 543             write (iout,*) 'Lorenz constants'
 544             do k=1,nlor(i,j)
 545               write (iout,'(3(1pe15.5))')
 546      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
 547             enddo
 548           enddo
 549         enddo
 550       endif
 551 C
 552 C 6/23/01 Read parameters for double torsionals
 553 C
 554       do i=1,ntortyp
 555         do j=1,ntortyp
 556           do k=1,ntortyp
 557             read (itordp,'(3a1)') t1,t2,t3
 558             if (t1.ne.onelett(i) .or. t2.ne.onelett(j)
 559      &        .or. t3.ne.onelett(k)) then
 560               write (iout,*) "Error in double torsional parameter file",
 561      &         i,j,k,t1,t2,t3
 562                stop "Error in double torsional parameter file"
 563             endif
 564             read (itordp,*) ntermd_1(i,j,k),ntermd_2(i,j,k)
 565             read (itordp,*) (v1c(1,l,i,j,k),l=1,ntermd_1(i,j,k))
 566             read (itordp,*) (v1s(1,l,i,j,k),l=1,ntermd_1(i,j,k))
 567             read (itordp,*) (v1c(2,l,i,j,k),l=1,ntermd_1(i,j,k))
 568             read (itordp,*) (v1s(2,l,i,j,k),l=1,ntermd_1(i,j,k))
 569             read (itordp,*) ((v2c(l,m,i,j,k),v2c(m,l,i,j,k),
 570      &       v2s(l,m,i,j,k),v2s(m,l,i,j,k),m=1,l-1),l=1,ntermd_2(i,j,k))
 571           enddo
 572         enddo
 573       enddo
 574       if (lprint) then
 575       write (iout,*)
 576       write (iout,*) 'Constants for double torsionals'
 577       do i=1,ntortyp
 578         do j=1,ntortyp
 579           do k=1,ntortyp
 580             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
 581      &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
 582             write (iout,*)
 583             write (iout,*) 'Single angles:'
 584             do l=1,ntermd_1(i,j,k)
 585               write (iout,'(i5,2f10.5,5x,2f10.5)') l,
 586      &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
 587      &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
 588             enddo
 589             write (iout,*)
 590             write (iout,*) 'Pairs of angles:'
 591             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 592             do l=1,ntermd_2(i,j,k)
 593               write (iout,'(i5,20f10.5)')
 594      &         l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 595             enddo
 596             write (iout,*)
 597             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
 598             do l=1,ntermd_2(i,j,k)
 599               write (iout,'(i5,20f10.5)')
 600      &         l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
 601             enddo
 602             write (iout,*)
 603           enddo
 604         enddo
 605       enddo
 606       endif
 607 #endif
 608 C Read of Side-chain backbone correlation parameters
 609 C Modified 11 May 2012 by Adasko
 610 CCC
 611 C
 612       read (isccor,*) nsccortyp
 613       read (isccor,*) (isccortyp(i),i=1,ntyp)
 614 c      write (iout,*) 'ntortyp',ntortyp
 615       maxinter=3
 616 cc maxinter is maximum interaction sites
 617       do l=1,maxinter
 618       do i=1,nsccortyp
 619         do j=1,nsccortyp
 620           read (isccor,*) nterm_sccor(i,j),nlor_sccor(i,j)
 621           v0ijsccor=0.0d0
 622           si=-1.0d0
 623
 624           do k=1,nterm_sccor(i,j)
 625             read (isccor,*) kk,v1sccor(k,l,i,j)
 626      &    ,v2sccor(k,l,i,j)
 627             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
 628             si=-si
 629           enddo
 630           do k=1,nlor_sccor(i,j)
 631             read (isccor,*) kk,vlor1sccor(k,i,j),
 632      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 633             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
 634      &(1+vlor3sccor(k,i,j)**2)
 635           enddo
 636           v0sccor(i,j)=v0ijsccor
 637         enddo
 638       enddo
 639       enddo
 640       close (isccor)
 641
 642       if (lprint) then
 643         write (iout,'(/a/)') 'Torsional constants:'
 644         do i=1,nsccortyp
 645           do j=1,nsccortyp
 646             write (iout,*) 'ityp',i,' jtyp',j
 647             write (iout,*) 'Fourier constants'
 648             do k=1,nterm_sccor(i,j)
 649            write (iout,'(2(1pe15.5))')v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 650             enddo
 651             write (iout,*) 'Lorenz constants'
 652             do k=1,nlor_sccor(i,j)
 653               write (iout,'(3(1pe15.5))')
 654      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 655             enddo
 656           enddo
 657         enddo
 658       endif
 659
 660 C
 661 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 662 C         interaction energy of the Gly, Ala, and Pro prototypes.
 663 C
 664       read (ifourier,*) nloctyp
 665       do i=1,nloctyp
 666         read (ifourier,*)
 667         read (ifourier,*) (b(ii,i),ii=1,13)
 668         if (lprint) then
 669         write (iout,*) 'Type',i
 670         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 671         endif
 672         B1(1,i)  = b(3,i)
 673         B1(2,i)  = b(5,i)
 674         B1tilde(1,i) = b(3,i)
 675         B1tilde(2,i) =-b(5,i)
 676         B2(1,i)  = b(2,i)
 677         B2(2,i)  = b(4,i)
 678         CC(1,1,i)= b(7,i)
 679         CC(2,2,i)=-b(7,i)
 680         CC(2,1,i)= b(9,i)
 681         CC(1,2,i)= b(9,i)
 682         Ctilde(1,1,i)=b(7,i)
 683         Ctilde(1,2,i)=b(9,i)
 684         Ctilde(2,1,i)=-b(9,i)
 685         Ctilde(2,2,i)=b(7,i)
 686         DD(1,1,i)= b(6,i)
 687         DD(2,2,i)=-b(6,i)
 688         DD(2,1,i)= b(8,i)
 689         DD(1,2,i)= b(8,i)
 690         Dtilde(1,1,i)=b(6,i)
 691         Dtilde(1,2,i)=b(8,i)
 692         Dtilde(2,1,i)=-b(8,i)
 693         Dtilde(2,2,i)=b(6,i)
 694         EE(1,1,i)= b(10,i)+b(11,i)
 695         EE(2,2,i)=-b(10,i)+b(11,i)
 696         EE(2,1,i)= b(12,i)-b(13,i)
 697         EE(1,2,i)= b(12,i)+b(13,i)
 698       enddo
 699       if (lprint) then
 700       do i=1,nloctyp
 701         write (iout,*) 'Type',i
 702         write (iout,*) 'B1'
 703 c        write (iout,'(f10.5)') B1(:,i)
 704         write(iout,*) B1(1,i),B1(2,i)
 705         write (iout,*) 'B2'
 706 c        write (iout,'(f10.5)') B2(:,i)
 707         write(iout,*) B2(1,i),B2(2,i)
 708         write (iout,*) 'CC'
 709         do j=1,2
 710           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 711         enddo
 712         write(iout,*) 'DD'
 713         do j=1,2
 714           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 715         enddo
 716         write(iout,*) 'EE'
 717         do j=1,2
 718           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 719         enddo
 720       enddo
 721       endif
 722 C
 723 C Read electrostatic-interaction parameters
 724 C
 725       if (lprint) then
 726         write (iout,'(/a)') 'Electrostatic interaction constants:'
 727         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 728      &            'IT','JT','APP','BPP','AEL6','AEL3'
 729       endif
 730       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
 731       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
 732       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
 733       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
 734       close (ielep)
 735       do i=1,2
 736         do j=1,2
 737         rri=rpp(i,j)**6
 738         app (i,j)=epp(i,j)*rri*rri
 739         bpp (i,j)=-2.0D0*epp(i,j)*rri
 740         ael6(i,j)=elpp6(i,j)*4.2D0**6
 741         ael3(i,j)=elpp3(i,j)*4.2D0**3
 742         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 743      &                    ael6(i,j),ael3(i,j)
 744         enddo
 745       enddo
 746 C
 747 C Read side-chain interaction parameters.
 748 C
 749       read (isidep,*) ipot,expon
 750       if (ipot.lt.1 .or. ipot.gt.5) then
 751         write (iout,'(2a)') 'Error while reading SC interaction',
 752      &               'potential file - unknown potential type.'
 753         stop
 754       endif
 755       expon2=expon/2
 756       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 757      & ', exponents are ',expon,2*expon
 758       goto (10,20,30,30,40) ipot
 759 C----------------------- LJ potential ---------------------------------
 760    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
 761       if (lprint) then
 762         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 763         write (iout,'(a/)') 'The epsilon array:'
 764         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 765         write (iout,'(/a)') 'One-body parameters:'
 766         write (iout,'(a,4x,a)') 'residue','sigma'
 767         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 768       endif
 769       goto 50
 770 C----------------------- LJK potential --------------------------------
 771    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 772      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 773       if (lprint) then
 774         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 775         write (iout,'(a/)') 'The epsilon array:'
 776         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 777         write (iout,'(/a)') 'One-body parameters:'
 778         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 779         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 780      &        i=1,ntyp)
 781       endif
 782       goto 50
 783 C---------------------- GB or BP potential -----------------------------
 784    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 785      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
 786      &  (alp(i),i=1,ntyp)
 787 C For the GB potential convert sigma'**2 into chi'
 788       if (ipot.eq.4) then
 789         do i=1,ntyp
 790           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
 791         enddo
 792       endif
 793       if (lprint) then
 794         write (iout,'(/a/)') 'Parameters of the BP potential:'
 795         write (iout,'(a/)') 'The epsilon array:'
 796         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 797         write (iout,'(/a)') 'One-body parameters:'
 798         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 799      &       '    chip  ','    alph  '
 800         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 801      &                     chip(i),alp(i),i=1,ntyp)
 802       endif
 803       goto 50
 804 C--------------------- GBV potential -----------------------------------
 805    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 806      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 807      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 808       if (lprint) then
 809         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 810         write (iout,'(a/)') 'The epsilon array:'
 811         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 812         write (iout,'(/a)') 'One-body parameters:'
 813         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 814      &      's||/s_|_^2','    chip  ','    alph  '
 815         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 816      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 817       endif
 818    50 continue
 819       close (isidep)
 820 C-----------------------------------------------------------------------
 821 C Calculate the "working" parameters of SC interactions.
 822       do i=2,ntyp
 823         do j=1,i-1
 824           eps(i,j)=eps(j,i)
 825         enddo
 826       enddo
 827       do i=1,ntyp
 828         do j=i,ntyp
 829           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 830           sigma(j,i)=sigma(i,j)
 831           rs0(i,j)=dwa16*sigma(i,j)
 832           rs0(j,i)=rs0(i,j)
 833         enddo
 834       enddo
 835       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 836      & 'Working parameters of the SC interactions:',
 837      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 838      & '  chi1   ','   chi2   '
 839       do i=1,ntyp
 840         do j=i,ntyp
 841           epsij=eps(i,j)
 842           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 843             rrij=sigma(i,j)
 844           else
 845             rrij=rr0(i)+rr0(j)
 846           endif
 847           r0(i,j)=rrij
 848           r0(j,i)=rrij
 849           rrij=rrij**expon
 850           epsij=eps(i,j)
 851           sigeps=dsign(1.0D0,epsij)
 852           epsij=dabs(epsij)
 853           aa(i,j)=epsij*rrij*rrij
 854           bb(i,j)=-sigeps*epsij*rrij
 855           aa(j,i)=aa(i,j)
 856           bb(j,i)=bb(i,j)
 857           if (ipot.gt.2) then
 858             sigt1sq=sigma0(i)**2
 859             sigt2sq=sigma0(j)**2
 860             sigii1=sigii(i)
 861             sigii2=sigii(j)
 862             ratsig1=sigt2sq/sigt1sq
 863             ratsig2=1.0D0/ratsig1
 864             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 865             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 866             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 867           else
 868             rsum_max=sigma(i,j)
 869           endif
 870 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 871             sigmaii(i,j)=rsum_max
 872             sigmaii(j,i)=rsum_max
 873 c         else
 874 c           sigmaii(i,j)=r0(i,j)
 875 c           sigmaii(j,i)=r0(i,j)
 876 c         endif
 877 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
 878           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
 879             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 880             augm(i,j)=epsij*r_augm**(2*expon)
 881 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 882             augm(j,i)=augm(i,j)
 883           else
 884             augm(i,j)=0.0D0
 885             augm(j,i)=0.0D0
 886           endif
 887           if (lprint) then
 888             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 889      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 890      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 891           endif
 892         enddo
 893       enddo
 894 C
 895 C Define the SC-p interaction constants and SS bond potentials
 896 C
 897       if (dyn_ss) then
 898         ss_depth=ebr/wsc-0.25*eps(1,1)
 899         Ht=Ht/wsc-0.25*eps(1,1)
 900         akcm=akcm*wstrain/wsc
 901         akth=akth*wstrain/wsc
 902         akct=akct*wstrain/wsc
 903         v1ss=v1ss*wstrain/wsc
 904         v2ss=v2ss*wstrain/wsc
 905         v3ss=v3ss*wstrain/wsc
 906       else
 907         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 908       endif
 909        write (iout,*) "Parameters of the SS-bond potential:"
 910        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 911      & " AKCT",akct
 912        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 913        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
 914        write (iout,*)" HT",Ht
 915
 916 #ifdef OLDSCP
 917       do i=1,20
 918 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
 919 C helix formation)
 920 c       aad(i,1)=0.3D0*4.0D0**12
 921 C Following line for constants currently implemented
 922 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
 923         aad(i,1)=1.5D0*4.0D0**12
 924 c       aad(i,1)=0.17D0*5.6D0**12
 925         aad(i,2)=aad(i,1)
 926 C "Soft" SC-p repulsion
 927         bad(i,1)=0.0D0
 928 C Following line for constants currently implemented
 929 c       aad(i,1)=0.3D0*4.0D0**6
 930 C "Hard" SC-p repulsion
 931         bad(i,1)=3.0D0*4.0D0**6
 932 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
 933         bad(i,2)=bad(i,1)
 934 c       aad(i,1)=0.0D0
 935 c       aad(i,2)=0.0D0
 936 c       bad(i,1)=1228.8D0
 937 c       bad(i,2)=1228.8D0
 938       enddo
 939 #else
 940 C
 941 C 8/9/01 Read the SC-p interaction constants from file
 942 C
 943       do i=1,ntyp
 944         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
 945       enddo
 946       do i=1,ntyp
 947         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
 948         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
 949         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
 950         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
 951       enddo
 952
 953       if (lprint) then
 954         write (iout,*) "Parameters of SC-p interactions:"
 955         do i=1,20
 956           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
 957      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
 958         enddo
 959       endif
 960 #endif
 961 C
 962 C Define the constants of the disulfide bridge
 963 C
 964 c      ebr=-5.50D0
 965 c
 966 c Old arbitrary potential - commented out.
 967 c
 968 c      dbr= 4.20D0
 969 c      fbr= 3.30D0
 970 c
 971 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
 972 c energy surface of diethyl disulfide.
 973 c A. Liwo and U. Kozlowska, 11/24/03
 974 c
 975 c      D0CM = 3.78d0
 976 c      AKCM = 15.1d0
 977 c      AKTH = 11.0d0
 978 c      AKCT = 12.0d0
 979 c      V1SS =-1.08d0
 980 c      V2SS = 7.61d0
 981 c      V3SS = 13.7d0
 982
 983 c      if (lprint) then
 984 c      write (iout,'(/a)') "Disulfide bridge parameters:"
 985 c      write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
 986 c      write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
 987 c      write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
 988 c      write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
 989 c     & ' v3ss:',v3ss
 990 c      endif
 991       return
 992       end