source/wham/src-restraints-DFA/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 i=1,nsccortyp
 651           do j=1,nsccortyp
 652             write (iout,*) 'ityp',i,' jtyp',j
 653             write (iout,*) 'Fourier constants'
 654             do k=1,nterm_sccor(i,j)
 655            write (iout,'(2(1pe15.5))')v1sccor(k,l,i,j),v2sccor(k,l,i,j)
 656             enddo
 657             write (iout,*) 'Lorenz constants'
 658             do k=1,nlor_sccor(i,j)
 659               write (iout,'(3(1pe15.5))')
 660      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
 661             enddo
 662           enddo
 663         enddo
 664       endif
 665
 666 C
 667 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 668 C         interaction energy of the Gly, Ala, and Pro prototypes.
 669 C
 670       read (ifourier,*) nloctyp
 671       do i=1,nloctyp
 672         read (ifourier,*)
 673         read (ifourier,*) (b(ii,i),ii=1,13)
 674         if (lprint) then
 675         write (iout,*) 'Type',i
 676         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 677         endif
 678         B1(1,i)  = b(3,i)
 679         B1(2,i)  = b(5,i)
 680         B1tilde(1,i) = b(3,i)
 681         B1tilde(2,i) =-b(5,i)
 682         B2(1,i)  = b(2,i)
 683         B2(2,i)  = b(4,i)
 684         CC(1,1,i)= b(7,i)
 685         CC(2,2,i)=-b(7,i)
 686         CC(2,1,i)= b(9,i)
 687         CC(1,2,i)= b(9,i)
 688         Ctilde(1,1,i)=b(7,i)
 689         Ctilde(1,2,i)=b(9,i)
 690         Ctilde(2,1,i)=-b(9,i)
 691         Ctilde(2,2,i)=b(7,i)
 692         DD(1,1,i)= b(6,i)
 693         DD(2,2,i)=-b(6,i)
 694         DD(2,1,i)= b(8,i)
 695         DD(1,2,i)= b(8,i)
 696         Dtilde(1,1,i)=b(6,i)
 697         Dtilde(1,2,i)=b(8,i)
 698         Dtilde(2,1,i)=-b(8,i)
 699         Dtilde(2,2,i)=b(6,i)
 700         EE(1,1,i)= b(10,i)+b(11,i)
 701         EE(2,2,i)=-b(10,i)+b(11,i)
 702         EE(2,1,i)= b(12,i)-b(13,i)
 703         EE(1,2,i)= b(12,i)+b(13,i)
 704       enddo
 705       if (lprint) then
 706       do i=1,nloctyp
 707         write (iout,*) 'Type',i
 708         write (iout,*) 'B1'
 709 c        write (iout,'(f10.5)') B1(:,i)
 710         write(iout,*) B1(1,i),B1(2,i)
 711         write (iout,*) 'B2'
 712 c        write (iout,'(f10.5)') B2(:,i)
 713         write(iout,*) B2(1,i),B2(2,i)
 714         write (iout,*) 'CC'
 715         do j=1,2
 716           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
 717         enddo
 718         write(iout,*) 'DD'
 719         do j=1,2
 720           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
 721         enddo
 722         write(iout,*) 'EE'
 723         do j=1,2
 724           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
 725         enddo
 726       enddo
 727       endif
 728 C
 729 C Read electrostatic-interaction parameters
 730 C
 731       if (lprint) then
 732         write (iout,'(/a)') 'Electrostatic interaction constants:'
 733         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
 734      &            'IT','JT','APP','BPP','AEL6','AEL3'
 735       endif
 736       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
 737       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
 738       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
 739       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
 740       close (ielep)
 741       do i=1,2
 742         do j=1,2
 743         rri=rpp(i,j)**6
 744         app (i,j)=epp(i,j)*rri*rri
 745         bpp (i,j)=-2.0D0*epp(i,j)*rri
 746         ael6(i,j)=elpp6(i,j)*4.2D0**6
 747         ael3(i,j)=elpp3(i,j)*4.2D0**3
 748         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
 749      &                    ael6(i,j),ael3(i,j)
 750         enddo
 751       enddo
 752 C
 753 C Read side-chain interaction parameters.
 754 C
 755       read (isidep,*) ipot,expon
 756       if (ipot.lt.1 .or. ipot.gt.5) then
 757         write (iout,'(2a)') 'Error while reading SC interaction',
 758      &               'potential file - unknown potential type.'
 759         stop
 760       endif
 761       expon2=expon/2
 762       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
 763      & ', exponents are ',expon,2*expon
 764       goto (10,20,30,30,40) ipot
 765 C----------------------- LJ potential ---------------------------------
 766    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
 767       if (lprint) then
 768         write (iout,'(/a/)') 'Parameters of the LJ potential:'
 769         write (iout,'(a/)') 'The epsilon array:'
 770         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 771         write (iout,'(/a)') 'One-body parameters:'
 772         write (iout,'(a,4x,a)') 'residue','sigma'
 773         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
 774       endif
 775       goto 50
 776 C----------------------- LJK potential --------------------------------
 777    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 778      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
 779       if (lprint) then
 780         write (iout,'(/a/)') 'Parameters of the LJK potential:'
 781         write (iout,'(a/)') 'The epsilon array:'
 782         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 783         write (iout,'(/a)') 'One-body parameters:'
 784         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
 785         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
 786      &        i=1,ntyp)
 787       endif
 788       goto 50
 789 C---------------------- GB or BP potential -----------------------------
 790    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 791      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
 792      &  (alp(i),i=1,ntyp)
 793 C For the GB potential convert sigma'**2 into chi'
 794       if (ipot.eq.4) then
 795         do i=1,ntyp
 796           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
 797         enddo
 798       endif
 799       if (lprint) then
 800         write (iout,'(/a/)') 'Parameters of the BP potential:'
 801         write (iout,'(a/)') 'The epsilon array:'
 802         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 803         write (iout,'(/a)') 'One-body parameters:'
 804         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
 805      &       '    chip  ','    alph  '
 806         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
 807      &                     chip(i),alp(i),i=1,ntyp)
 808       endif
 809       goto 50
 810 C--------------------- GBV potential -----------------------------------
 811    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
 812      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
 813      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
 814       if (lprint) then
 815         write (iout,'(/a/)') 'Parameters of the GBV potential:'
 816         write (iout,'(a/)') 'The epsilon array:'
 817         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
 818         write (iout,'(/a)') 'One-body parameters:'
 819         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
 820      &      's||/s_|_^2','    chip  ','    alph  '
 821         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
 822      &           sigii(i),chip(i),alp(i),i=1,ntyp)
 823       endif
 824    50 continue
 825       close (isidep)
 826 C-----------------------------------------------------------------------
 827 C Calculate the "working" parameters of SC interactions.
 828       do i=2,ntyp
 829         do j=1,i-1
 830           eps(i,j)=eps(j,i)
 831         enddo
 832       enddo
 833       do i=1,ntyp
 834         do j=i,ntyp
 835           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
 836           sigma(j,i)=sigma(i,j)
 837           rs0(i,j)=dwa16*sigma(i,j)
 838           rs0(j,i)=rs0(i,j)
 839         enddo
 840       enddo
 841       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
 842      & 'Working parameters of the SC interactions:',
 843      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
 844      & '  chi1   ','   chi2   '
 845       do i=1,ntyp
 846         do j=i,ntyp
 847           epsij=eps(i,j)
 848           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 849             rrij=sigma(i,j)
 850           else
 851             rrij=rr0(i)+rr0(j)
 852           endif
 853           r0(i,j)=rrij
 854           r0(j,i)=rrij
 855           rrij=rrij**expon
 856           epsij=eps(i,j)
 857           sigeps=dsign(1.0D0,epsij)
 858           epsij=dabs(epsij)
 859           aa(i,j)=epsij*rrij*rrij
 860           bb(i,j)=-sigeps*epsij*rrij
 861           aa(j,i)=aa(i,j)
 862           bb(j,i)=bb(i,j)
 863           if (ipot.gt.2) then
 864             sigt1sq=sigma0(i)**2
 865             sigt2sq=sigma0(j)**2
 866             sigii1=sigii(i)
 867             sigii2=sigii(j)
 868             ratsig1=sigt2sq/sigt1sq
 869             ratsig2=1.0D0/ratsig1
 870             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
 871             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
 872             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
 873           else
 874             rsum_max=sigma(i,j)
 875           endif
 876 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
 877             sigmaii(i,j)=rsum_max
 878             sigmaii(j,i)=rsum_max
 879 c         else
 880 c           sigmaii(i,j)=r0(i,j)
 881 c           sigmaii(j,i)=r0(i,j)
 882 c         endif
 883 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
 884           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
 885             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
 886             augm(i,j)=epsij*r_augm**(2*expon)
 887 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
 888             augm(j,i)=augm(i,j)
 889           else
 890             augm(i,j)=0.0D0
 891             augm(j,i)=0.0D0
 892           endif
 893           if (lprint) then
 894             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
 895      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
 896      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
 897           endif
 898         enddo
 899       enddo
 900 C
 901 C Define the SC-p interaction constants and SS bond potentials
 902 C
 903       if (dyn_ss) then
 904         ss_depth=ebr/wsc-0.25*eps(1,1)
 905         Ht=Ht/wsc-0.25*eps(1,1)
 906         akcm=akcm*wstrain/wsc
 907         akth=akth*wstrain/wsc
 908         akct=akct*wstrain/wsc
 909         v1ss=v1ss*wstrain/wsc
 910         v2ss=v2ss*wstrain/wsc
 911         v3ss=v3ss*wstrain/wsc
 912       else
 913         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 914       endif
 915        write (iout,*) "Parameters of the SS-bond potential:"
 916        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 917      & " AKCT",akct
 918        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 919        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
 920        write (iout,*)" HT",Ht
 921
 922 #ifdef OLDSCP
 923       do i=1,20
 924 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
 925 C helix formation)
 926 c       aad(i,1)=0.3D0*4.0D0**12
 927 C Following line for constants currently implemented
 928 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
 929         aad(i,1)=1.5D0*4.0D0**12
 930 c       aad(i,1)=0.17D0*5.6D0**12
 931         aad(i,2)=aad(i,1)
 932 C "Soft" SC-p repulsion
 933         bad(i,1)=0.0D0
 934 C Following line for constants currently implemented
 935 c       aad(i,1)=0.3D0*4.0D0**6
 936 C "Hard" SC-p repulsion
 937         bad(i,1)=3.0D0*4.0D0**6
 938 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
 939         bad(i,2)=bad(i,1)
 940 c       aad(i,1)=0.0D0
 941 c       aad(i,2)=0.0D0
 942 c       bad(i,1)=1228.8D0
 943 c       bad(i,2)=1228.8D0
 944       enddo
 945 #else
 946 C
 947 C 8/9/01 Read the SC-p interaction constants from file
 948 C
 949       do i=1,ntyp
 950         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
 951       enddo
 952       do i=1,ntyp
 953         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
 954         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
 955         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
 956         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
 957       enddo
 958
 959       if (lprint) then
 960         write (iout,*) "Parameters of SC-p interactions:"
 961         do i=1,20
 962           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
 963      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
 964         enddo
 965       endif
 966 #endif
 967 C
 968 C Define the constants of the disulfide bridge
 969 C
 970 c      ebr=-5.50D0
 971 c
 972 c Old arbitrary potential - commented out.
 973 c
 974 c      dbr= 4.20D0
 975 c      fbr= 3.30D0
 976 c
 977 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
 978 c energy surface of diethyl disulfide.
 979 c A. Liwo and U. Kozlowska, 11/24/03
 980 c
 981 c      D0CM = 3.78d0
 982 c      AKCM = 15.1d0
 983 c      AKTH = 11.0d0
 984 c      AKCT = 12.0d0
 985 c      V1SS =-1.08d0
 986 c      V2SS = 7.61d0
 987 c      V3SS = 13.7d0
 988
 989 c      if (lprint) then
 990 c      write (iout,'(/a)') "Disulfide bridge parameters:"
 991 c      write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
 992 c      write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
 993 c      write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
 994 c      write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
 995 c     & ' v3ss:',v3ss
 996 c      endif
 997       return
 998       end