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