source/wham/src-M-SAXS/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       include 'COMMON.SHIELD'
  25       include 'COMMON.CONTROL'
  26       include 'COMMON.LANGEVIN'
  27       character*1 t1,t2,t3
  28       character*1 onelett(4) /"G","A","P","D"/
  29       character*1 toronelet(-2:2) /"p","a","G","A","P"/
  30       logical lprint
  31       dimension blower(3,3,maxlob)
  32       character*800 controlcard
  33       character*256 bondname_t,thetname_t,rotname_t,torname_t,
  34      & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
  35      & sccorname_t
  36       integer ilen
  37       external ilen
  38       character*16 key
  39       integer iparm
  40       character*6 res1
  41       character*3 lancuch,ucase
  42 C      write (iout,*) "KURWA"
  43 C
  44 C Body
  45 C
  46       call getenv("PRINT_PARM",lancuch)
  47       lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
  48       write (iout,*) "lprint ",lprint
  49 C Set LPRINT=.TRUE. for debugging
  50       dwa16=2.0d0**(1.0d0/6.0d0)
  51       itypro=20
  52 C Assign virtual-bond length
  53       vbl=3.8D0
  54       vblinv=1.0D0/vbl
  55       vblinv2=vblinv*vblinv
  56       call card_concat(controlcard,.true.)
  57       wname(4)="WCORRH"
  58       do i=1,n_ene
  59         key = wname(i)(:ilen(wname(i)))
  60         call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
  61       enddo
  62
  63       write (iout,*) "iparm",iparm," myparm",myparm
  64 c If reading not own parameters, skip assignment
  65       call reada(controlcard,"D0CM",d0cm,3.78d0)
  66       call reada(controlcard,"AKCM",akcm,15.1d0)
  67       call reada(controlcard,"AKTH",akth,11.0d0)
  68       call reada(controlcard,"AKCT",akct,12.0d0)
  69       call reada(controlcard,"V1SS",v1ss,-1.08d0)
  70       call reada(controlcard,"V2SS",v2ss,7.61d0)
  71       call reada(controlcard,"V3SS",v3ss,13.7d0)
  72       call reada(controlcard,"EBR",ebr,-5.50D0)
  73       call reada(controlcard,"DTRISS",dtriss,1.0D0)
  74       call reada(controlcard,"ATRISS",atriss,0.3D0)
  75       call reada(controlcard,"BTRISS",btriss,0.02D0)
  76       call reada(controlcard,"CTRISS",ctriss,1.0D0)
  77       dyn_ss=(index(controlcard,'DYN_SS').gt.0)
  78       write(iout,*) "ATRISS",atriss
  79       write(iout,*) "BTRISS",btriss
  80       write(iout,*) "CTRISS",ctriss
  81       write(iout,*) "DTRISS",dtriss
  82
  83 C      do i=1,maxres
  84 C        dyn_ss_mask(i)=.false.
  85 C      enddo
  86 C      ebr=-12.0D0
  87 c
  88 c Old arbitrary potential - commented out.
  89 c
  90 c      dbr= 4.20D0
  91 c      fbr= 3.30D0
  92 c
  93 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
  94 c energy surface of diethyl disulfide.
  95 c A. Liwo and U. Kozlowska, 11/24/03
  96 c
  97       D0CM = 3.78d0
  98       AKCM = 15.1d0
  99       AKTH = 11.0d0
 100       AKCT = 12.0d0
 101       V1SS =-1.08d0
 102       V2SS = 7.61d0
 103       V3SS = 13.7d0
 104
 105       do i=1,maxres-1
 106         do j=i+1,maxres
 107           dyn_ssbond_ij(i,j)=1.0d300
 108         enddo
 109       enddo
 110       call reada(controlcard,"HT",Ht,0.0D0)
 111 C      if (dyn_ss) then
 112 C        ss_depth=ebr/wsc-0.25*eps(1,1)
 113 C        write(iout,*) HT,wsc,eps(1,1),'KURWA'
 114 C        Ht=Ht/wsc-0.25*eps(1,1)
 115
 116 C        akcm=akcm*whpb/wsc
 117 C        akth=akth*whpb/wsc
 118 C        akct=akct*whpb/wsc
 119 C        v1ss=v1ss*whpb/wsc
 120 C        v2ss=v2ss*whpb/wsc
 121 C        v3ss=v3ss*whpb/wsc
 122 C      else
 123 C        ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
 124 C      endif
 125
 126       if (iparm.eq.myparm .or. .not.separate_parset) then
 127
 128 c
 129 c Setup weights for UNRES
 130 c
 131       wsc=ww(1)
 132       wscp=ww(2)
 133       welec=ww(3)
 134       wcorr=ww(4)
 135       wcorr5=ww(5)
 136       wcorr6=ww(6)
 137       wel_loc=ww(7)
 138       wturn3=ww(8)
 139       wturn4=ww(9)
 140       wturn6=ww(10)
 141       wang=ww(11)
 142       wscloc=ww(12)
 143       wtor=ww(13)
 144       wtor_d=ww(14)
 145       wvdwpp=ww(16)
 146       wbond=ww(18)
 147       wsccor=ww(19)
 148       whpb=ww(15)
 149       wstrain=ww(15)
 150       wliptran=ww(22)
 151       wshield=ww(25)
 152       wsaxs=ww(26)
 153       write (iout,*) "PARMREAD: wsaxs",wsaxs
 154       endif
 155
 156       call card_concat(controlcard,.false.)
 157
 158 c Return if not own parameters
 159
 160       if (iparm.ne.myparm .and. separate_parset) return
 161
 162       call reads(controlcard,"BONDPAR",bondname_t,bondname)
 163       open (ibond,file=bondname_t,status='old')
 164       rewind(ibond)
 165       call reads(controlcard,"THETPAR",thetname_t,thetname)
 166       open (ithep,file=thetname_t,status='old')
 167       rewind(ithep)
 168       call reads(controlcard,"ROTPAR",rotname_t,rotname)
 169       open (irotam,file=rotname_t,status='old')
 170       rewind(irotam)
 171       call reads(controlcard,"TORPAR",torname_t,torname)
 172       open (itorp,file=torname_t,status='old')
 173       rewind(itorp)
 174       call reads(controlcard,"TORDPAR",tordname_t,tordname)
 175       write (iout,*) "tor_mode",tor_mode
 176       call flush(iout)
 177       if (tor_mode.eq.0)
 178      & open (itordp,file=tordname_t,status='old')
 179       rewind(itordp)
 180       call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
 181       open (isccor,file=sccorname_t,status='old')
 182       rewind(isccor)
 183       call reads(controlcard,"FOURIER",fouriername_t,fouriername)
 184       open (ifourier,file=fouriername_t,status='old')
 185       rewind(ifourier)
 186       call reads(controlcard,"ELEPAR",elename_t,elename)
 187       open (ielep,file=elename_t,status='old')
 188       rewind(ielep)
 189       call reads(controlcard,"SIDEPAR",sidename_t,sidename)
 190       open (isidep,file=sidename_t,status='old')
 191       rewind(isidep)
 192       call reads(controlcard,"SCPPAR",scpname_t,scpname)
 193       open (iscpp,file=scpname_t,status='old')
 194       rewind(iscpp)
 195       write (iout,*) "Parameter set:",iparm
 196       write (iout,*) "Energy-term weights:"
 197       do i=1,n_ene
 198         write (iout,'(a16,f10.5)') wname(i),ww(i)
 199       enddo
 200       write (iout,*) "Sidechain potential file        : ",
 201      &  sidename_t(:ilen(sidename_t))
 202 #ifndef OLDSCP
 203       write (iout,*) "SCp potential file              : ",
 204      &  scpname_t(:ilen(scpname_t))
 205 #endif
 206       write (iout,*) "Electrostatic potential file    : ",
 207      &  elename_t(:ilen(elename_t))
 208       write (iout,*) "Cumulant coefficient file       : ",
 209      &  fouriername_t(:ilen(fouriername_t))
 210       write (iout,*) "Torsional parameter file        : ",
 211      &  torname_t(:ilen(torname_t))
 212       write (iout,*) "Double torsional parameter file : ",
 213      &  tordname_t(:ilen(tordname_t))
 214       write (iout,*) "Backbone-rotamer parameter file : ",
 215      &  sccorname(:ilen(sccorname))
 216       write (iout,*) "Bond & inertia constant file    : ",
 217      &  bondname_t(:ilen(bondname_t))
 218       write (iout,*) "Bending parameter file          : ",
 219      &  thetname_t(:ilen(thetname_t))
 220       write (iout,*) "Rotamer parameter file          : ",
 221      &  rotname_t(:ilen(rotname_t))
 222
 223 c
 224 c Read the virtual-bond parameters, masses, and moments of inertia
 225 c and Stokes' radii of the peptide group and side chains
 226 c
 227 #ifdef CRYST_BOND
 228       read (ibond,*,end=121,err=121) vbldp0,vbldpdum,akp,mp,ip,pstok
 229       do i=1,ntyp
 230         nbondterm(i)=1
 231         read (ibond,*,end=121,err=121) vbldsc0(1,i),aksc(1,i),
 232      &    msc(i),isc(i),restok(i)
 233         dsc(i) = vbldsc0(1,i)
 234         if (i.eq.10) then
 235           dsc_inv(i)=0.0D0
 236         else
 237           dsc_inv(i)=1.0D0/dsc(i)
 238         endif
 239       enddo
 240 #else
 241       read (ibond,*,end=121,err=121) ijunk,vbldp0,vbldpdum,akp,rjunk,
 242      &   mp,ip,pstok
 243       do i=1,ntyp
 244         read (ibond,*,end=121,err=121) nbondterm(i),(vbldsc0(j,i),
 245      &   aksc(j,i),abond0(j,i),j=1,nbondterm(i)),msc(i),isc(i),restok(i)
 246         dsc(i) = vbldsc0(1,i)
 247         if (i.eq.10) then
 248           dsc_inv(i)=0.0D0
 249         else
 250           dsc_inv(i)=1.0D0/dsc(i)
 251         endif
 252       enddo
 253 #endif
 254       if (lprint) then
 255         write(iout,'(/a/)')"Force constants virtual bonds:"
 256         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
 257      &   'inertia','Pstok'
 258         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
 259         do i=1,ntyp
 260           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
 261      &      vbldsc0(1,i),aksc(1,i),abond0(1,i)
 262           do j=2,nbondterm(i)
 263             write (iout,'(13x,3f10.5)')
 264      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
 265           enddo
 266         enddo
 267       endif
 268        write (iout,*) "iliptranpar",iliptranpar
 269        write (iout,*) "liptranname ",liptranname
 270        read(iliptranpar,*,end=1161,err=1161) pepliptran
 271        write (iout,*) "pepliptran",pepliptran
 272        do i=1,ntyp
 273        read(iliptranpar,*,end=1161,err=1161) liptranene(i)
 274        write (iout,*) i,liptranene(i)
 275        enddo
 276        rewind iliptranpar
 277 #ifdef CRYST_THETA
 278 C
 279 C Read the parameters of the probability distribution/energy expression
 280 C of the virtual-bond valence angles theta
 281 C
 282       do i=1,ntyp
 283         read (ithep,*,end=111,err=111) a0thet(i),(athet(j,i,1,1),j=1,2),
 284      &    (bthet(j,i,1,1),j=1,2)
 285         read (ithep,*,end=111,err=111) (polthet(j,i),j=0,3)
 286         read (ithep,*,end=111,err=111) (gthet(j,i),j=1,3)
 287         read (ithep,*,end=111,err=111) theta0(i),sig0(i),sigc0(i)
 288         sigc0(i)=sigc0(i)**2
 289       enddo
 290       do i=1,ntyp
 291       athet(1,i,1,-1)=athet(1,i,1,1)
 292       athet(2,i,1,-1)=athet(2,i,1,1)
 293       bthet(1,i,1,-1)=-bthet(1,i,1,1)
 294       bthet(2,i,1,-1)=-bthet(2,i,1,1)
 295       athet(1,i,-1,1)=-athet(1,i,1,1)
 296       athet(2,i,-1,1)=-athet(2,i,1,1)
 297       bthet(1,i,-1,1)=bthet(1,i,1,1)
 298       bthet(2,i,-1,1)=bthet(2,i,1,1)
 299       enddo
 300       do i=-ntyp,-1
 301       a0thet(i)=a0thet(-i)
 302       athet(1,i,-1,-1)=athet(1,-i,1,1)
 303       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 304       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 305       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 306       athet(1,i,-1,1)=athet(1,-i,1,1)
 307       athet(2,i,-1,1)=-athet(2,-i,1,1)
 308       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 309       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 310       athet(1,i,1,-1)=-athet(1,-i,1,1)
 311       athet(2,i,1,-1)=athet(2,-i,1,1)
 312       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 313       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 314       theta0(i)=theta0(-i)
 315       sig0(i)=sig0(-i)
 316       sigc0(i)=sigc0(-i)
 317        do j=0,3
 318         polthet(j,i)=polthet(j,-i)
 319        enddo
 320        do j=1,3
 321          gthet(j,i)=gthet(j,-i)
 322        enddo
 323       enddo
 324       close (ithep)
 325       if (lprint) then
 326 c       write (iout,'(a)')
 327 c    &    'Parameters of the virtual-bond valence angles:'
 328 c       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
 329 c    & '    ATHETA0   ','         A1   ','        A2    ',
 330 c    & '        B1    ','         B2   '
 331 c       do i=1,ntyp
 332 c         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 333 c    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
 334 c       enddo
 335 c       write (iout,'(/a/9x,5a/79(1h-))')
 336 c    & 'Parameters of the expression for sigma(theta_c):',
 337 c    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
 338 c    & '     ALPH3    ','    SIGMA0C   '
 339 c       do i=1,ntyp
 340 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 341 c    &      (polthet(j,i),j=0,3),sigc0(i)
 342 c       enddo
 343 c       write (iout,'(/a/9x,5a/79(1h-))')
 344 c    & 'Parameters of the second gaussian:',
 345 c    & '    THETA0    ','     SIGMA0   ','        G1    ',
 346 c    & '        G2    ','         G3   '
 347 c       do i=1,ntyp
 348 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
 349 c    &       sig0(i),(gthet(j,i),j=1,3)
 350 c       enddo
 351         write (iout,'(a)')
 352      &    'Parameters of the virtual-bond valence angles:'
 353         write (iout,'(/a/9x,5a/79(1h-))')
 354      & 'Coefficients of expansion',
 355      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
 356      & '   b1*10^1    ','    b2*10^1   '
 357         do i=1,ntyp
 358           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
 359      &        a0thet(i),(100*athet(j,i,1,1),j=1,2),
 360      &        (10*bthet(j,i,1,1),j=1,2)
 361         enddo
 362         write (iout,'(/a/9x,5a/79(1h-))')
 363      & 'Parameters of the expression for sigma(theta_c):',
 364      & ' alpha0       ','  alph1       ',' alph2        ',
 365      & ' alhp3        ','   sigma0c    '
 366         do i=1,ntyp
 367           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
 368      &      (polthet(j,i),j=0,3),sigc0(i)
 369         enddo
 370         write (iout,'(/a/9x,5a/79(1h-))')
 371      & 'Parameters of the second gaussian:',
 372      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
 373      & '        G2    ','   G3*10^1    '
 374         do i=1,ntyp
 375           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
 376      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 377         enddo
 378       endif
 379 #else
 380       IF (tor_mode.eq.0) THEN
 381 C
 382 C Read the parameters of Utheta determined from ab initio surfaces
 383 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 384 C
 385       read (ithep,*,end=111,err=111) nthetyp,ntheterm,ntheterm2,
 386      &  ntheterm3,nsingle,ndouble
 387       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 388       read (ithep,*,end=111,err=111) (ithetyp(i),i=1,ntyp1)
 389       do i=-ntyp1,-1
 390         ithetyp(i)=-ithetyp(-i)
 391       enddo
 392       do iblock=1,2
 393       do i=-maxthetyp,maxthetyp
 394         do j=-maxthetyp,maxthetyp
 395           do k=-maxthetyp,maxthetyp
 396             aa0thet(i,j,k,iblock)=0.0d0
 397             do l=1,ntheterm
 398               aathet(l,i,j,k,iblock)=0.0d0
 399             enddo
 400             do l=1,ntheterm2
 401               do m=1,nsingle
 402                 bbthet(m,l,i,j,k,iblock)=0.0d0
 403                 ccthet(m,l,i,j,k,iblock)=0.0d0
 404                 ddthet(m,l,i,j,k,iblock)=0.0d0
 405                 eethet(m,l,i,j,k,iblock)=0.0d0
 406               enddo
 407             enddo
 408             do l=1,ntheterm3
 409               do m=1,ndouble
 410                 do mm=1,ndouble
 411                  ffthet(mm,m,l,i,j,k,iblock)=0.0d0
 412                  ggthet(mm,m,l,i,j,k,iblock)=0.0d0
 413                 enddo
 414               enddo
 415             enddo
 416           enddo
 417         enddo
 418       enddo
 419       enddo
 420 C      write (iout,*) "KURWA1"
 421       do iblock=1,2
 422       do i=0,nthetyp
 423         do j=-nthetyp,nthetyp
 424           do k=-nthetyp,nthetyp
 425             read (ithep,'(6a)',end=111,err=111) res1
 426             write(iout,*) res1,i,j,k
 427             read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
 428             read (ithep,*,end=111,err=111)(aathet(l,i,j,k,iblock),
 429      &        l=1,ntheterm)
 430             read (ithep,*,end=111,err=111)
 431      &       ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 432      &        (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 433      &        (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),
 434      &        (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle)
 435      &        ,ll=1,ntheterm2)
 436             read (ithep,*,end=111,err=111)
 437      &      (((ffthet(llll,lll,ll,i,j,k,iblock),
 438      &      ffthet(lll,llll,ll,i,j,k,iblock),
 439      &         ggthet(llll,lll,ll,i,j,k,iblock)
 440      &        ,ggthet(lll,llll,ll,i,j,k,iblock),
 441      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 442           enddo
 443         enddo
 444       enddo
 445 C       write(iout,*) "KURWA1.1"
 446 C
 447 C For dummy ends assign glycine-type coefficients of theta-only terms; the
 448 C coefficients of theta-and-gamma-dependent terms are zero.
 449 C
 450       do i=1,nthetyp
 451         do j=1,nthetyp
 452           do l=1,ntheterm
 453             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
 454             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
 455           enddo
 456           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
 457           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
 458         enddo
 459         do l=1,ntheterm
 460           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 461         enddo
 462         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 463       enddo
 464       enddo
 465 C       write(iout,*) "KURWA1.5"
 466 C Substitution for D aminoacids from symmetry.
 467       do iblock=1,2
 468       do i=-nthetyp,0
 469         do j=-nthetyp,nthetyp
 470           do k=-nthetyp,nthetyp
 471            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
 472            do l=1,ntheterm
 473            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
 474            enddo
 475            do ll=1,ntheterm2
 476             do lll=1,nsingle
 477             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
 478             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
 479             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
 480             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
 481             enddo
 482           enddo
 483           do ll=1,ntheterm3
 484            do lll=2,ndouble
 485             do llll=1,lll-1
 486             ffthet(llll,lll,ll,i,j,k,iblock)=
 487      &      ffthet(llll,lll,ll,-i,-j,-k,iblock)
 488             ffthet(lll,llll,ll,i,j,k,iblock)=
 489      &      ffthet(lll,llll,ll,-i,-j,-k,iblock)
 490             ggthet(llll,lll,ll,i,j,k,iblock)=
 491      &      -ggthet(llll,lll,ll,-i,-j,-k,iblock)
 492             ggthet(lll,llll,ll,i,j,k,iblock)=
 493      &      -ggthet(lll,llll,ll,-i,-j,-k,iblock)
 494             enddo !ll
 495            enddo  !lll
 496           enddo   !llll
 497          enddo    !k
 498         enddo     !j
 499        enddo      !i
 500       enddo       !iblock
 501
 502 C
 503 C Control printout of the coefficients of virtual-bond-angle potentials
 504 C
 505       if (lprint) then
 506         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 507         do iblock=1,2
 508         do i=1,nthetyp+1
 509           do j=1,nthetyp+1
 510             do k=1,nthetyp+1
 511               write (iout,'(//4a)')
 512      &         'Type ',onelett(i),onelett(j),onelett(k)
 513               write (iout,'(//a,10x,a)') " l","a[l]"
 514               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
 515               write (iout,'(i2,1pe15.5)')
 516      &           (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
 517             do l=1,ntheterm2
 518               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
 519      &          "b",l,"c",l,"d",l,"e",l
 520               do m=1,nsingle
 521                 write (iout,'(i2,4(1pe15.5))') m,
 522      &          bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),
 523      &          ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
 524               enddo
 525             enddo
 526             do l=1,ntheterm3
 527               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
 528      &          "f+",l,"f-",l,"g+",l,"g-",l
 529               do m=2,ndouble
 530                 do n=1,m-1
 531                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
 532      &              ffthet(n,m,l,i,j,k,iblock),
 533      &              ffthet(m,n,l,i,j,k,iblock),
 534      &              ggthet(n,m,l,i,j,k,iblock),
 535      &              ggthet(m,n,l,i,j,k,iblock)
 536                 enddo
 537               enddo
 538             enddo
 539           enddo
 540         enddo
 541       enddo
 542       enddo
 543       call flush(iout)
 544       endif
 545
 546       ELSE
 547
 548 C here will be the apropriate recalibrating for D-aminoacid
 549       read (ithep,*,end=111,err=111) nthetyp
 550       do i=-nthetyp+1,nthetyp-1
 551         read (ithep,*,end=111,err=111) nbend_kcc_Tb(i)
 552         do j=0,nbend_kcc_Tb(i)
 553           read (ithep,*,end=111,err=111) ijunk,v1bend_chyb(j,i)
 554         enddo
 555       enddo
 556       if (lprint) then
 557         write (iout,'(a)')
 558      &    "Parameters of the valence-only potentials"
 559         do i=-nthetyp+1,nthetyp-1
 560         write (iout,'(2a)') "Type ",toronelet(i)
 561         do k=0,nbend_kcc_Tb(i)
 562           write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
 563         enddo
 564         enddo
 565       endif
 566
 567       ENDIF ! TOR_MODE
 568
 569       close(ithep)
 570 #endif
 571 C      write(iout,*) 'KURWA2'
 572 #ifdef CRYST_SC
 573 C
 574 C Read the parameters of the probability distribution/energy expression
 575 C of the side chains.
 576 C
 577       do i=1,ntyp
 578 cc      write (iout,*) "tu dochodze",i
 579         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
 580         if (i.eq.10) then
 581           dsc_inv(i)=0.0D0
 582         else
 583           dsc_inv(i)=1.0D0/dsc(i)
 584         endif
 585         if (i.ne.10) then
 586         do j=1,nlob(i)
 587           do k=1,3
 588             do l=1,3
 589               blower(l,k,j)=0.0D0
 590             enddo
 591           enddo
 592         enddo
 593         bsc(1,i)=0.0D0
 594         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
 595      &    ((blower(k,l,1),l=1,k),k=1,3)
 596         censc(1,1,-i)=censc(1,1,i)
 597         censc(2,1,-i)=censc(2,1,i)
 598         censc(3,1,-i)=-censc(3,1,i)
 599         do j=2,nlob(i)
 600           read (irotam,*,end=112,err=112) bsc(j,i)
 601           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
 602      &                                 ((blower(k,l,j),l=1,k),k=1,3)
 603         censc(1,j,-i)=censc(1,j,i)
 604         censc(2,j,-i)=censc(2,j,i)
 605         censc(3,j,-i)=-censc(3,j,i)
 606 C BSC is amplitude of Gaussian
 607         enddo
 608         do j=1,nlob(i)
 609           do k=1,3
 610             do l=1,k
 611               akl=0.0D0
 612               do m=1,3
 613                 akl=akl+blower(k,m,j)*blower(l,m,j)
 614               enddo
 615               gaussc(k,l,j,i)=akl
 616               gaussc(l,k,j,i)=akl
 617              if (((k.eq.3).and.(l.ne.3))
 618      &        .or.((l.eq.3).and.(k.ne.3))) then
 619                 gaussc(k,l,j,-i)=-akl
 620                 gaussc(l,k,j,-i)=-akl
 621               else
 622                 gaussc(k,l,j,-i)=akl
 623                 gaussc(l,k,j,-i)=akl
 624               endif
 625             enddo
 626           enddo
 627         enddo
 628         endif
 629       enddo
 630       close (irotam)
 631       if (lprint) then
 632         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 633         do i=1,ntyp
 634           nlobi=nlob(i)
 635           if (nlobi.gt.0) then
 636           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
 637      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 638 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 639 c          write (iout,'(a,f10.4,4(16x,f10.4))')
 640 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
 641 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
 642            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
 643      &                             'log h',(bsc(j,i),j=1,nlobi)
 644            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
 645      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 646 c          write (iout,'(a)')
 647 c         do j=1,nlobi
 648 c           ind=0
 649 c           do k=1,3
 650 c             do l=1,k
 651 c              ind=ind+1
 652 c              blower(k,l,j)=gaussc(ind,j,i)
 653 c             enddo
 654 c           enddo
 655 c         enddo
 656           do k=1,3
 657             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
 658      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 659           enddo
 660           endif
 661         enddo
 662       endif
 663 #else
 664 C
 665 C Read scrot parameters for potentials determined from all-atom AM1 calculations
 666 C added by Urszula Kozlowska 07/11/2007
 667 C
 668       do i=1,ntyp
 669         read (irotam,*,end=112,err=112)
 670        if (i.eq.10) then
 671          read (irotam,*,end=112,err=112)
 672        else
 673          do j=1,65
 674            read(irotam,*,end=112,err=112) sc_parmin(j,i)
 675          enddo
 676        endif
 677       enddo
 678 #endif
 679       close(irotam)
 680 C
 681 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 682 C         interaction energy of the Gly, Ala, and Pro prototypes.
 683 C
 684       read (ifourier,*,end=115,err=115) nloctyp
 685       SPLIT_FOURIERTOR = nloctyp.lt.0
 686       nloctyp = iabs(nloctyp)
 687 #ifdef NEWCORR
 688       read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
 689       read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
 690       itype2loc(ntyp1)=nloctyp
 691       iloctyp(nloctyp)=ntyp1
 692       do i=1,ntyp1
 693         itype2loc(-i)=-itype2loc(i)
 694       enddo
 695 #else
 696       iloctyp(0)=10
 697       iloctyp(1)=9
 698       iloctyp(2)=20
 699       iloctyp(3)=ntyp1
 700 #endif
 701       do i=1,nloctyp
 702         iloctyp(-i)=-iloctyp(i)
 703       enddo
 704 c      write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
 705 c      write (iout,*) "nloctyp",nloctyp,
 706 c     &  " iloctyp",(iloctyp(i),i=0,nloctyp)
 707 #ifdef NEWCORR
 708       do i=0,nloctyp-1
 709 c             write (iout,*) "NEWCORR",i
 710         read (ifourier,*,end=115,err=115)
 711         do ii=1,3
 712           do j=1,2
 713             read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
 714           enddo
 715         enddo
 716 c             write (iout,*) "NEWCORR BNEW1"
 717 c             write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
 718         do ii=1,3
 719           do j=1,2
 720             read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
 721           enddo
 722         enddo
 723 c             write (iout,*) "NEWCORR BNEW2"
 724 c             write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
 725         do kk=1,3
 726           read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
 727           read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
 728         enddo
 729 c             write (iout,*) "NEWCORR CCNEW"
 730 c             write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
 731         do kk=1,3
 732           read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
 733           read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
 734         enddo
 735 c             write (iout,*) "NEWCORR DDNEW"
 736 c             write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
 737         do ii=1,2
 738           do jj=1,2
 739             do kk=1,2
 740               read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
 741             enddo
 742           enddo
 743         enddo
 744 c             write (iout,*) "NEWCORR EENEW1"
 745 c             write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
 746         do ii=1,3
 747           read (ifourier,*,end=115,err=115) e0new(ii,i)
 748         enddo
 749 c             write (iout,*) (e0new(ii,i),ii=1,3)
 750       enddo
 751 c             write (iout,*) "NEWCORR EENEW"
 752       do i=0,nloctyp-1
 753         do ii=1,3
 754           ccnew(ii,1,i)=ccnew(ii,1,i)/2
 755           ccnew(ii,2,i)=ccnew(ii,2,i)/2
 756           ddnew(ii,1,i)=ddnew(ii,1,i)/2
 757           ddnew(ii,2,i)=ddnew(ii,2,i)/2
 758         enddo
 759       enddo
 760       do i=1,nloctyp-1
 761         do ii=1,3
 762           bnew1(ii,1,-i)= bnew1(ii,1,i)
 763           bnew1(ii,2,-i)=-bnew1(ii,2,i)
 764           bnew2(ii,1,-i)= bnew2(ii,1,i)
 765           bnew2(ii,2,-i)=-bnew2(ii,2,i)
 766         enddo
 767         do ii=1,3
 768 c          ccnew(ii,1,i)=ccnew(ii,1,i)/2
 769 c          ccnew(ii,2,i)=ccnew(ii,2,i)/2
 770 c          ddnew(ii,1,i)=ddnew(ii,1,i)/2
 771 c          ddnew(ii,2,i)=ddnew(ii,2,i)/2
 772           ccnew(ii,1,-i)=ccnew(ii,1,i)
 773           ccnew(ii,2,-i)=-ccnew(ii,2,i)
 774           ddnew(ii,1,-i)=ddnew(ii,1,i)
 775           ddnew(ii,2,-i)=-ddnew(ii,2,i)
 776         enddo
 777         e0new(1,-i)= e0new(1,i)
 778         e0new(2,-i)=-e0new(2,i)
 779         e0new(3,-i)=-e0new(3,i)
 780         do kk=1,2
 781           eenew(kk,1,1,-i)= eenew(kk,1,1,i)
 782           eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
 783           eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
 784           eenew(kk,2,2,-i)= eenew(kk,2,2,i)
 785         enddo
 786       enddo
 787       if (lprint) then
 788         write (iout,'(a)') "Coefficients of the multibody terms"
 789         do i=-nloctyp+1,nloctyp-1
 790           write (iout,*) "Type: ",onelet(iloctyp(i))
 791           write (iout,*) "Coefficients of the expansion of B1"
 792           do j=1,2
 793             write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
 794           enddo
 795           write (iout,*) "Coefficients of the expansion of B2"
 796           do j=1,2
 797             write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
 798           enddo
 799           write (iout,*) "Coefficients of the expansion of C"
 800           write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
 801           write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
 802           write (iout,*) "Coefficients of the expansion of D"
 803           write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
 804           write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
 805           write (iout,*) "Coefficients of the expansion of E"
 806           write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
 807           do j=1,2
 808             do k=1,2
 809               write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
 810             enddo
 811           enddo
 812         enddo
 813       endif
 814       IF (SPLIT_FOURIERTOR) THEN
 815       do i=0,nloctyp-1
 816 c             write (iout,*) "NEWCORR TOR",i
 817         read (ifourier,*,end=115,err=115)
 818         do ii=1,3
 819           do j=1,2
 820             read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
 821           enddo
 822         enddo
 823 c             write (iout,*) "NEWCORR BNEW1 TOR"
 824 c             write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
 825         do ii=1,3
 826           do j=1,2
 827             read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
 828           enddo
 829         enddo
 830 c             write (iout,*) "NEWCORR BNEW2 TOR"
 831 c             write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
 832         do kk=1,3
 833           read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
 834           read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
 835         enddo
 836 c             write (iout,*) "NEWCORR CCNEW TOR"
 837 c             write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
 838         do kk=1,3
 839           read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
 840           read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
 841         enddo
 842 c             write (iout,*) "NEWCORR DDNEW TOR"
 843 c             write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
 844         do ii=1,2
 845           do jj=1,2
 846             do kk=1,2
 847               read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
 848             enddo
 849           enddo
 850         enddo
 851 c         write (iout,*) "NEWCORR EENEW1 TOR"
 852 c         write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
 853         do ii=1,3
 854           read (ifourier,*,end=115,err=115) e0newtor(ii,i)
 855         enddo
 856 c             write (iout,*) (e0newtor(ii,i),ii=1,3)
 857       enddo
 858 c             write (iout,*) "NEWCORR EENEW TOR"
 859       do i=0,nloctyp-1
 860         do ii=1,3
 861           ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
 862           ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
 863           ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
 864           ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
 865         enddo
 866       enddo
 867       do i=1,nloctyp-1
 868         do ii=1,3
 869           bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
 870           bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
 871           bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
 872           bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
 873         enddo
 874         do ii=1,3
 875           ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
 876           ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
 877           ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
 878           ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
 879         enddo
 880         e0newtor(1,-i)= e0newtor(1,i)
 881         e0newtor(2,-i)=-e0newtor(2,i)
 882         e0newtor(3,-i)=-e0newtor(3,i)
 883         do kk=1,2
 884           eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
 885           eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
 886           eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
 887           eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
 888         enddo
 889       enddo
 890       if (lprint) then
 891         write (iout,'(a)')
 892      &    "Single-body coefficients of the torsional potentials"
 893         do i=-nloctyp+1,nloctyp-1
 894           write (iout,*) "Type: ",onelet(iloctyp(i))
 895           write (iout,*) "Coefficients of the expansion of B1tor"
 896           do j=1,2
 897             write (iout,'(3hB1(,i1,1h),3f10.5)')
 898      &        j,(bnew1tor(k,j,i),k=1,3)
 899           enddo
 900           write (iout,*) "Coefficients of the expansion of B2tor"
 901           do j=1,2
 902             write (iout,'(3hB2(,i1,1h),3f10.5)')
 903      &        j,(bnew2tor(k,j,i),k=1,3)
 904           enddo
 905           write (iout,*) "Coefficients of the expansion of Ctor"
 906           write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
 907           write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
 908           write (iout,*) "Coefficients of the expansion of Dtor"
 909           write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
 910           write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
 911           write (iout,*) "Coefficients of the expansion of Etor"
 912           write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
 913           do j=1,2
 914             do k=1,2
 915               write (iout,'(1hE,2i1,2f10.5)')
 916      &          j,k,(eenewtor(l,j,k,i),l=1,2)
 917             enddo
 918           enddo
 919         enddo
 920       endif
 921       ELSE
 922       do i=-nloctyp+1,nloctyp-1
 923         do ii=1,3
 924           do j=1,2
 925             bnew1tor(ii,j,i)=bnew1(ii,j,i)
 926           enddo
 927         enddo
 928         do ii=1,3
 929           do j=1,2
 930             bnew2tor(ii,j,i)=bnew2(ii,j,i)
 931           enddo
 932         enddo
 933         do ii=1,3
 934           ccnewtor(ii,1,i)=ccnew(ii,1,i)
 935           ccnewtor(ii,2,i)=ccnew(ii,2,i)
 936           ddnewtor(ii,1,i)=ddnew(ii,1,i)
 937           ddnewtor(ii,2,i)=ddnew(ii,2,i)
 938         enddo
 939       enddo
 940       ENDIF !SPLIT_FOURIER_TOR
 941 #else
 942       if (lprint)
 943      &  write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
 944       do i=0,nloctyp-1
 945         read (ifourier,*,end=115,err=115)
 946         read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
 947         if (lprint) then
 948         write (iout,*) 'Type ',onelet(iloctyp(i))
 949         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
 950         endif
 951         if (i.gt.0) then
 952         b(2,-i)= b(2,i)
 953         b(3,-i)= b(3,i)
 954         b(4,-i)=-b(4,i)
 955         b(5,-i)=-b(5,i)
 956         endif
 957 c        B1(1,i)  = b(3)
 958 c        B1(2,i)  = b(5)
 959 c        B1(1,-i) = b(3)
 960 c        B1(2,-i) = -b(5)
 961 c        b1(1,i)=0.0d0
 962 c        b1(2,i)=0.0d0
 963 c        B1tilde(1,i) = b(3)
 964 c        B1tilde(2,i) =-b(5)
 965 c        B1tilde(1,-i) =-b(3)
 966 c        B1tilde(2,-i) =b(5)
 967 c        b1tilde(1,i)=0.0d0
 968 c        b1tilde(2,i)=0.0d0
 969 c        B2(1,i)  = b(2)
 970 c        B2(2,i)  = b(4)
 971 c        B2(1,-i)  =b(2)
 972 c        B2(2,-i)  =-b(4)
 973 cc        B1tilde(1,i) = b(3,i)
 974 cc        B1tilde(2,i) =-b(5,i)
 975 C        B1tilde(1,-i) =-b(3,i)
 976 C        B1tilde(2,-i) =b(5,i)
 977 cc        b1tilde(1,i)=0.0d0
 978 cc        b1tilde(2,i)=0.0d0
 979 cc        B2(1,i)  = b(2,i)
 980 cc        B2(2,i)  = b(4,i)
 981 C        B2(1,-i)  =b(2,i)
 982 C        B2(2,-i)  =-b(4,i)
 983
 984 c        b2(1,i)=0.0d0
 985 c        b2(2,i)=0.0d0
 986         CCold(1,1,i)= b(7,i)
 987         CCold(2,2,i)=-b(7,i)
 988         CCold(2,1,i)= b(9,i)
 989         CCold(1,2,i)= b(9,i)
 990         CCold(1,1,-i)= b(7,i)
 991         CCold(2,2,-i)=-b(7,i)
 992         CCold(2,1,-i)=-b(9,i)
 993         CCold(1,2,-i)=-b(9,i)
 994 c        CC(1,1,i)=0.0d0
 995 c        CC(2,2,i)=0.0d0
 996 c        CC(2,1,i)=0.0d0
 997 c        CC(1,2,i)=0.0d0
 998 c        Ctilde(1,1,i)= CCold(1,1,i)
 999 c        Ctilde(1,2,i)= CCold(1,2,i)
1000 c        Ctilde(2,1,i)=-CCold(2,1,i)
1001 c        Ctilde(2,2,i)=-CCold(2,2,i)
1002
1003 c        Ctilde(1,1,i)=0.0d0
1004 c        Ctilde(1,2,i)=0.0d0
1005 c        Ctilde(2,1,i)=0.0d0
1006 c        Ctilde(2,2,i)=0.0d0
1007         DDold(1,1,i)= b(6,i)
1008         DDold(2,2,i)=-b(6,i)
1009         DDold(2,1,i)= b(8,i)
1010         DDold(1,2,i)= b(8,i)
1011         DDold(1,1,-i)= b(6,i)
1012         DDold(2,2,-i)=-b(6,i)
1013         DDold(2,1,-i)=-b(8,i)
1014         DDold(1,2,-i)=-b(8,i)
1015 c        DD(1,1,i)=0.0d0
1016 c        DD(2,2,i)=0.0d0
1017 c        DD(2,1,i)=0.0d0
1018 c        DD(1,2,i)=0.0d0
1019 c        Dtilde(1,1,i)= DD(1,1,i)
1020 c        Dtilde(1,2,i)= DD(1,2,i)
1021 c        Dtilde(2,1,i)=-DD(2,1,i)
1022 c        Dtilde(2,2,i)=-DD(2,2,i)
1023
1024 c        Dtilde(1,1,i)=0.0d0
1025 c        Dtilde(1,2,i)=0.0d0
1026 c        Dtilde(2,1,i)=0.0d0
1027 c        Dtilde(2,2,i)=0.0d0
1028         EEold(1,1,i)= b(10,i)+b(11,i)
1029         EEold(2,2,i)=-b(10,i)+b(11,i)
1030         EEold(2,1,i)= b(12,i)-b(13,i)
1031         EEold(1,2,i)= b(12,i)+b(13,i)
1032         EEold(1,1,-i)= b(10,i)+b(11,i)
1033         EEold(2,2,-i)=-b(10,i)+b(11,i)
1034         EEold(2,1,-i)=-b(12,i)+b(13,i)
1035         EEold(1,2,-i)=-b(12,i)-b(13,i)
1036         write(iout,*) "TU DOCHODZE"
1037         print *,"JESTEM"
1038 c        ee(1,1,i)=1.0d0
1039 c        ee(2,2,i)=1.0d0
1040 c        ee(2,1,i)=0.0d0
1041 c        ee(1,2,i)=0.0d0
1042 c        ee(2,1,i)=ee(1,2,i)
1043       enddo
1044       if (lprint) then
1045       write (iout,*)
1046       write (iout,*)
1047      &"Coefficients of the cumulants (independent of valence angles)"
1048       do i=-nloctyp+1,nloctyp-1
1049         write (iout,*) 'Type ',onelet(iloctyp(i))
1050         write (iout,*) 'B1'
1051         write(iout,'(2f10.5)') B(3,i),B(5,i)
1052         write (iout,*) 'B2'
1053         write(iout,'(2f10.5)') B(2,i),B(4,i)
1054         write (iout,*) 'CC'
1055         do j=1,2
1056           write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1057         enddo
1058         write(iout,*) 'DD'
1059         do j=1,2
1060           write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1061         enddo
1062         write(iout,*) 'EE'
1063         do j=1,2
1064           write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1065         enddo
1066       enddo
1067       endif
1068 #endif
1069 C      write (iout,*) 'KURWAKURWA'
1070 #ifdef CRYST_TOR
1071 C
1072 C Read torsional parameters in old format
1073 C
1074       read (itorp,*,end=113,err=113) ntortyp,nterm_old
1075       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1076       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1077       do i=1,ntortyp
1078         do j=1,ntortyp
1079           read (itorp,'(a)',end=113,err=113)
1080           do k=1,nterm_old
1081             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1082           enddo
1083         enddo
1084       enddo
1085       close (itorp)
1086       if (lprint) then
1087         write (iout,'(/a/)') 'Torsional constants:'
1088         do i=1,ntortyp
1089           do j=1,ntortyp
1090             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1091             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1092           enddo
1093         enddo
1094       endif
1095 #else
1096 C
1097 C Read torsional parameters
1098 C
1099       IF (TOR_MODE.eq.0) THEN
1100
1101       read (itorp,*,end=113,err=113) ntortyp
1102       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1103       do i=1,ntyp1
1104         itype2loc(i)=itortyp(i)
1105       enddo
1106       write (iout,*) 'ntortyp',ntortyp
1107       do i=1,ntyp1
1108         itype2loc(-i)=-itype2loc(i)
1109       enddo
1110       itortyp(ntyp1)=ntortyp
1111       do iblock=1,2
1112       do i=-ntyp,-1
1113        itortyp(i)=-itortyp(-i)
1114       enddo
1115 c      write (iout,*) 'ntortyp',ntortyp
1116       do i=0,ntortyp-1
1117         do j=-ntortyp+1,ntortyp-1
1118           read (itorp,*,end=113,err=113) nterm(i,j,iblock),
1119      &          nlor(i,j,iblock)
1120           nterm(-i,-j,iblock)=nterm(i,j,iblock)
1121           nlor(-i,-j,iblock)=nlor(i,j,iblock)
1122           v0ij=0.0d0
1123           si=-1.0d0
1124           do k=1,nterm(i,j,iblock)
1125             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),
1126      &      v2(k,i,j,iblock)
1127             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1128             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1129             v0ij=v0ij+si*v1(k,i,j,iblock)
1130             si=-si
1131          enddo
1132           do k=1,nlor(i,j,iblock)
1133             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
1134      &        vlor2(k,i,j),vlor3(k,i,j)
1135             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1136           enddo
1137           v0(i,j,iblock)=v0ij
1138           v0(-i,-j,iblock)=v0ij
1139         enddo
1140       enddo
1141       enddo
1142       close (itorp)
1143       if (lprint) then
1144         write (iout,'(/a/)') 'Torsional constants:'
1145         do i=1,ntortyp
1146           do j=1,ntortyp
1147             write (iout,*) 'ityp',i,' jtyp',j
1148             write (iout,*) 'Fourier constants'
1149             do k=1,nterm(i,j,iblock)
1150               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),
1151      &        v2(k,i,j,iblock)
1152             enddo
1153             write (iout,*) 'Lorenz constants'
1154             do k=1,nlor(i,j,iblock)
1155               write (iout,'(3(1pe15.5))')
1156      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1157             enddo
1158           enddo
1159         enddo
1160       endif
1161 C
1162 C 6/23/01 Read parameters for double torsionals
1163 C
1164       do iblock=1,2
1165       do i=0,ntortyp-1
1166         do j=-ntortyp+1,ntortyp-1
1167           do k=-ntortyp+1,ntortyp-1
1168             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1169 c              write (iout,*) "OK onelett",
1170 c     &         i,j,k,t1,t2,t3
1171
1172             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j)
1173      &        .or. t3.ne.toronelet(k)) then
1174               write (iout,*) "Error in double torsional parameter file",
1175      &         i,j,k,t1,t2,t3
1176 #ifdef MPI
1177               call MPI_Finalize(Ierror)
1178 #endif
1179                stop "Error in double torsional parameter file"
1180             endif
1181           read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),
1182      &         ntermd_2(i,j,k,iblock)
1183             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1184             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1185             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,
1186      &         ntermd_1(i,j,k,iblock))
1187             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,
1188      &         ntermd_1(i,j,k,iblock))
1189             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,
1190      &         ntermd_1(i,j,k,iblock))
1191             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,
1192      &         ntermd_1(i,j,k,iblock))
1193 C Martix of D parameters for one dimesional foureir series
1194             do l=1,ntermd_1(i,j,k,iblock)
1195              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1196              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1197              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1198              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1199 c            write(iout,*) "whcodze" ,
1200 c     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1201             enddo
1202             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),
1203      &         v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),
1204      &         v2s(m,l,i,j,k,iblock),
1205      &         m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1206 C Martix of D parameters for two dimesional fourier series
1207             do l=1,ntermd_2(i,j,k,iblock)
1208              do m=1,l-1
1209              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1210              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1211              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1212              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1213              enddo!m
1214             enddo!l
1215           enddo!k
1216         enddo!j
1217       enddo!i
1218       enddo!iblock
1219       if (lprint) then
1220       write (iout,*)
1221       write (iout,*) 'Constants for double torsionals'
1222       do iblock=1,2
1223       do i=0,ntortyp-1
1224         do j=-ntortyp+1,ntortyp-1
1225           do k=-ntortyp+1,ntortyp-1
1226             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
1227      &        ' nsingle',ntermd_1(i,j,k,iblock),
1228      &        ' ndouble',ntermd_2(i,j,k,iblock)
1229             write (iout,*)
1230             write (iout,*) 'Single angles:'
1231             do l=1,ntermd_1(i,j,k,iblock)
1232               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,
1233      &           v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),
1234      &           v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),
1235      &           v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1236             enddo
1237             write (iout,*)
1238             write (iout,*) 'Pairs of angles:'
1239             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1240             do l=1,ntermd_2(i,j,k,iblock)
1241               write (iout,'(i5,20f10.5)')
1242      &         l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1243             enddo
1244             write (iout,*)
1245            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1246             do l=1,ntermd_2(i,j,k,iblock)
1247               write (iout,'(i5,20f10.5)')
1248      &         l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),
1249      &         (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1250             enddo
1251             write (iout,*)
1252           enddo
1253         enddo
1254       enddo
1255       enddo
1256       endif
1257
1258       ELSE IF (TOR_MODE.eq.1) THEN
1259
1260 C read valence-torsional parameters
1261       read (itorp,*,end=113,err=113) ntortyp
1262       nkcctyp=ntortyp
1263       write (iout,*) "Valence-torsional parameters read in ntortyp",
1264      &   ntortyp
1265       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1266       write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1267       do i=-ntyp,-1
1268         itortyp(i)=-itortyp(-i)
1269       enddo
1270       do i=-ntortyp+1,ntortyp-1
1271         do j=-ntortyp+1,ntortyp-1
1272 C first we read the cos and sin gamma parameters
1273           read (itorp,'(13x,a)',end=113,err=113) string
1274           write (iout,*) i,j,string
1275           read (itorp,*,end=113,err=113)
1276      &    nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1277 C           read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1278           do k=1,nterm_kcc(j,i)
1279             do l=1,nterm_kcc_Tb(j,i)
1280               do ll=1,nterm_kcc_Tb(j,i)
1281               read (itorp,*,end=113,err=113) ii,jj,kk,
1282      &          v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1283               enddo
1284             enddo
1285           enddo
1286         enddo
1287       enddo
1288
1289       ELSE
1290 c AL 4/8/16: Calculate coefficients from one-body parameters
1291       ntortyp=nloctyp
1292       do i=-ntyp1,ntyp1
1293         itortyp(i)=itype2loc(i)
1294       enddo
1295       write (iout,*)
1296      &"Val-tor parameters calculated from cumulant coefficients ntortyp"
1297      & ,ntortyp
1298       do i=-ntortyp+1,ntortyp-1
1299         do j=-ntortyp+1,ntortyp-1
1300           nterm_kcc(j,i)=2
1301           nterm_kcc_Tb(j,i)=3
1302           do k=1,nterm_kcc_Tb(j,i)
1303             do l=1,nterm_kcc_Tb(j,i)
1304               v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
1305      &                         +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1306               v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
1307      &                         +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1308             enddo
1309           enddo
1310           do k=1,nterm_kcc_Tb(j,i)
1311             do l=1,nterm_kcc_Tb(j,i)
1312 #ifdef CORRCD
1313               v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1314      &                         -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1315               v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1316      &                         +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1317 #else
1318               v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1319      &                         -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1320               v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1321      &                         +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1322 #endif
1323             enddo
1324           enddo
1325 cf(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
1326         enddo
1327       enddo
1328
1329       ENDIF ! TOR_MODE
1330
1331       if (tor_mode.gt.0 .and. lprint) then
1332 c Print valence-torsional parameters
1333         write (iout,'(a)')
1334      &    "Parameters of the valence-torsional potentials"
1335         do i=-ntortyp+1,ntortyp-1
1336         do j=-ntortyp+1,ntortyp-1
1337         write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1338         write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1339         do k=1,nterm_kcc(j,i)
1340           do l=1,nterm_kcc_Tb(j,i)
1341             do ll=1,nterm_kcc_Tb(j,i)
1342                write (iout,'(3i5,2f15.4)')
1343      &            k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1344             enddo
1345           enddo
1346         enddo
1347         enddo
1348         enddo
1349       endif
1350
1351 #endif
1352 C Read of Side-chain backbone correlation parameters
1353 C Modified 11 May 2012 by Adasko
1354 CCC
1355 C
1356       read (isccor,*,end=119,err=119) nsccortyp
1357       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1358       do i=-ntyp,-1
1359         isccortyp(i)=-isccortyp(-i)
1360       enddo
1361       iscprol=isccortyp(20)
1362 c      write (iout,*) 'ntortyp',ntortyp
1363       maxinter=3
1364 cc maxinter is maximum interaction sites
1365       do l=1,maxinter
1366       do i=1,nsccortyp
1367         do j=1,nsccortyp
1368           read (isccor,*,end=119,err=119)
1369      &nterm_sccor(i,j),nlor_sccor(i,j)
1370 c          write (iout,*) nterm_sccor(i,j)
1371           v0ijsccor=0.0d0
1372           v0ijsccor1=0.0d0
1373           v0ijsccor2=0.0d0
1374           v0ijsccor3=0.0d0
1375           si=-1.0d0
1376           nterm_sccor(-i,j)=nterm_sccor(i,j)
1377           nterm_sccor(-i,-j)=nterm_sccor(i,j)
1378           nterm_sccor(i,-j)=nterm_sccor(i,j)
1379 c          write (iout,*) nterm_sccor(i,j),nterm_sccor(-i,j),
1380 c     &    nterm_sccor(-i,-j),nterm_sccor(i,-j)
1381           do k=1,nterm_sccor(i,j)
1382             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
1383      &    ,v2sccor(k,l,i,j)
1384             if (j.eq.iscprol) then
1385              if (i.eq.isccortyp(10)) then
1386              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1387              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1388              else
1389              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
1390      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1391              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
1392      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1393              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1394              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1395              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1396              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1397              endif
1398             else
1399              if (i.eq.isccortyp(10)) then
1400              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1401              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1402              else
1403                if (j.eq.isccortyp(10)) then
1404              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1405              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1406                else
1407              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1408              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1409              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1410              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1411              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1412              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1413                 endif
1414                endif
1415             endif
1416             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1417             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1418             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1419             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1420             si=-si
1421            enddo
1422           do k=1,nlor_sccor(i,j)
1423             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),
1424      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1425             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
1426      &(1+vlor3sccor(k,i,j)**2)
1427           enddo
1428           v0sccor(l,i,j)=v0ijsccor
1429           v0sccor(l,-i,j)=v0ijsccor1
1430           v0sccor(l,i,-j)=v0ijsccor2
1431           v0sccor(l,-i,-j)=v0ijsccor3
1432           enddo
1433         enddo
1434       enddo
1435       close (isccor)
1436       if (lprint) then
1437         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1438         do l=1,maxinter
1439         do i=1,nsccortyp
1440           do j=1,nsccortyp
1441             write (iout,*) 'ityp',i,' jtyp',j
1442             write (iout,*) 'Fourier constants'
1443             do k=1,nterm_sccor(i,j)
1444               write (iout,'(2(1pe15.5))')
1445      & v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1446             enddo
1447             write (iout,*) 'Lorenz constants'
1448             do k=1,nlor_sccor(i,j)
1449               write (iout,'(3(1pe15.5))')
1450      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1451             enddo
1452           enddo
1453         enddo
1454         enddo
1455       endif
1456 C
1457 C Read electrostatic-interaction parameters
1458 C
1459       if (lprint) then
1460         write (iout,'(/a)') 'Electrostatic interaction constants:'
1461         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1462      &            'IT','JT','APP','BPP','AEL6','AEL3'
1463       endif
1464       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1465       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1466       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1467       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1468       close (ielep)
1469       do i=1,2
1470         do j=1,2
1471         rri=rpp(i,j)**6
1472         app (i,j)=epp(i,j)*rri*rri
1473         bpp (i,j)=-2.0D0*epp(i,j)*rri
1474         ael6(i,j)=elpp6(i,j)*4.2D0**6
1475         ael3(i,j)=elpp3(i,j)*4.2D0**3
1476         lprint=.true.
1477         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1478      &                    ael6(i,j),ael3(i,j)
1479          lprint=.false.
1480         enddo
1481       enddo
1482 C
1483 C Read side-chain interaction parameters.
1484 C
1485       read (isidep,*,end=117,err=117) ipot,expon
1486       if (ipot.lt.1 .or. ipot.gt.5) then
1487         write (iout,'(2a)') 'Error while reading SC interaction',
1488      &               'potential file - unknown potential type.'
1489         stop
1490       endif
1491       expon2=expon/2
1492       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
1493      & ', exponents are ',expon,2*expon
1494       goto (10,20,30,30,40) ipot
1495 C----------------------- LJ potential ---------------------------------
1496    10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1497      &   (sigma0(i),i=1,ntyp)
1498       if (lprint) then
1499         write (iout,'(/a/)') 'Parameters of the LJ potential:'
1500         write (iout,'(a/)') 'The epsilon array:'
1501         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1502         write (iout,'(/a)') 'One-body parameters:'
1503         write (iout,'(a,4x,a)') 'residue','sigma'
1504         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1505       endif
1506       goto 50
1507 C----------------------- LJK potential --------------------------------
1508    20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1509      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1510       if (lprint) then
1511         write (iout,'(/a/)') 'Parameters of the LJK potential:'
1512         write (iout,'(a/)') 'The epsilon array:'
1513         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1514         write (iout,'(/a)') 'One-body parameters:'
1515         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1516         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
1517      &        i=1,ntyp)
1518       endif
1519       goto 50
1520 C---------------------- GB or BP potential -----------------------------
1521    30 do i=1,ntyp
1522        read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1523       enddo
1524       read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1525       read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1526       read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1527       read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1528       do i=1,ntyp
1529        read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
1530 C       write(iout,*) "WARNING!!",i,ntyp
1531        write(iout,*) "epslip", i, (epslip(i,j),j=i,ntyp)
1532 C       do j=1,ntyp
1533 C       epslip(i,j)=epslip(i,j)+0.05d0
1534 C       enddo
1535       enddo
1536 C For the GB potential convert sigma'**2 into chi'
1537       if (ipot.eq.4) then
1538         do i=1,ntyp
1539           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1540         enddo
1541       endif
1542       if (lprint) then
1543         write (iout,'(/a/)') 'Parameters of the BP potential:'
1544         write (iout,'(a/)') 'The epsilon array:'
1545         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1546         write (iout,'(/a)') 'One-body parameters:'
1547         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
1548      &       '    chip  ','    alph  '
1549         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
1550      &                     chip(i),alp(i),i=1,ntyp)
1551       endif
1552       goto 50
1553 C--------------------- GBV potential -----------------------------------
1554    40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),
1555      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
1556      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1557       if (lprint) then
1558         write (iout,'(/a/)') 'Parameters of the GBV potential:'
1559         write (iout,'(a/)') 'The epsilon array:'
1560         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1561         write (iout,'(/a)') 'One-body parameters:'
1562         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
1563      &      's||/s_|_^2','    chip  ','    alph  '
1564         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
1565      &           sigii(i),chip(i),alp(i),i=1,ntyp)
1566       endif
1567    50 continue
1568       close (isidep)
1569 C-----------------------------------------------------------------------
1570 C Calculate the "working" parameters of SC interactions.
1571       do i=2,ntyp
1572         do j=1,i-1
1573           eps(i,j)=eps(j,i)
1574           epslip(i,j)=epslip(j,i)
1575         enddo
1576       enddo
1577       do i=1,ntyp
1578         do j=i,ntyp
1579           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1580           sigma(j,i)=sigma(i,j)
1581           rs0(i,j)=dwa16*sigma(i,j)
1582           rs0(j,i)=rs0(i,j)
1583         enddo
1584       enddo
1585       if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1586      & 'Working parameters of the SC interactions:',
1587      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
1588      & '  chi1   ','   chi2   '
1589       do i=1,ntyp
1590         do j=i,ntyp
1591           epsij=eps(i,j)
1592           epsijlip=epslip(i,j)
1593           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1594             rrij=sigma(i,j)
1595           else
1596             rrij=rr0(i)+rr0(j)
1597           endif
1598           r0(i,j)=rrij
1599           r0(j,i)=rrij
1600           rrij=rrij**expon
1601           epsij=eps(i,j)
1602           sigeps=dsign(1.0D0,epsij)
1603           epsij=dabs(epsij)
1604           aa_aq(i,j)=epsij*rrij*rrij
1605           bb_aq(i,j)=-sigeps*epsij*rrij
1606           aa_aq(j,i)=aa_aq(i,j)
1607           bb_aq(j,i)=bb_aq(i,j)
1608           sigeps=dsign(1.0D0,epsijlip)
1609           epsijlip=dabs(epsijlip)
1610           aa_lip(i,j)=epsijlip*rrij*rrij
1611           bb_lip(i,j)=-sigeps*epsijlip*rrij
1612           aa_lip(j,i)=aa_lip(i,j)
1613           bb_lip(j,i)=bb_lip(i,j)
1614           if (ipot.gt.2) then
1615             sigt1sq=sigma0(i)**2
1616             sigt2sq=sigma0(j)**2
1617             sigii1=sigii(i)
1618             sigii2=sigii(j)
1619             ratsig1=sigt2sq/sigt1sq
1620             ratsig2=1.0D0/ratsig1
1621             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1622             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1623             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1624           else
1625             rsum_max=sigma(i,j)
1626           endif
1627 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1628             sigmaii(i,j)=rsum_max
1629             sigmaii(j,i)=rsum_max
1630 c         else
1631 c           sigmaii(i,j)=r0(i,j)
1632 c           sigmaii(j,i)=r0(i,j)
1633 c         endif
1634 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1635           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1636             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1637             augm(i,j)=epsij*r_augm**(2*expon)
1638 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1639             augm(j,i)=augm(i,j)
1640           else
1641             augm(i,j)=0.0D0
1642             augm(j,i)=0.0D0
1643           endif
1644           if (lprint) then
1645             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1646      &      restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
1647      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1648           endif
1649         enddo
1650       enddo
1651 C
1652 C Define the SC-p interaction constants
1653 C
1654 #ifdef OLDSCP
1655       do i=1,20
1656 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1657 C helix formation)
1658 c       aad(i,1)=0.3D0*4.0D0**12
1659 C Following line for constants currently implemented
1660 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
1661         aad(i,1)=1.5D0*4.0D0**12
1662 c       aad(i,1)=0.17D0*5.6D0**12
1663         aad(i,2)=aad(i,1)
1664 C "Soft" SC-p repulsion
1665         bad(i,1)=0.0D0
1666 C Following line for constants currently implemented
1667 c       aad(i,1)=0.3D0*4.0D0**6
1668 C "Hard" SC-p repulsion
1669         bad(i,1)=3.0D0*4.0D0**6
1670 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1671         bad(i,2)=bad(i,1)
1672 c       aad(i,1)=0.0D0
1673 c       aad(i,2)=0.0D0
1674 c       bad(i,1)=1228.8D0
1675 c       bad(i,2)=1228.8D0
1676       enddo
1677 #else
1678 C
1679 C 8/9/01 Read the SC-p interaction constants from file
1680 C
1681       do i=1,ntyp
1682         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
1683       enddo
1684       do i=1,ntyp
1685         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1686         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1687         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1688         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1689       enddo
1690
1691       if (lprint) then
1692         write (iout,*) "Parameters of SC-p interactions:"
1693         do i=1,20
1694           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1695      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1696         enddo
1697       endif
1698 #endif
1699 C
1700 C Define the constants of the disulfide bridge
1701 C
1702 C      ebr=-12.0D0
1703 c
1704 c Old arbitrary potential - commented out.
1705 c
1706 c      dbr= 4.20D0
1707 c      fbr= 3.30D0
1708 c
1709 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1710 c energy surface of diethyl disulfide.
1711 c A. Liwo and U. Kozlowska, 11/24/03
1712 c
1713 C      D0CM = 3.78d0
1714 C      AKCM = 15.1d0
1715 C      AKTH = 11.0d0
1716 C      AKCT = 12.0d0
1717 C      V1SS =-1.08d0
1718 C      V2SS = 7.61d0
1719 C      V3SS = 13.7d0
1720       write (iout,*) dyn_ss,'dyndyn'
1721       if (dyn_ss) then
1722         ss_depth=ebr/wsc-0.25*eps(1,1)
1723 C        write(iout,*) akcm,whpb,wsc,'KURWA'
1724         Ht=Ht/wsc-0.25*eps(1,1)
1725
1726         akcm=akcm*whpb/wsc
1727         akth=akth*whpb/wsc
1728         akct=akct*whpb/wsc
1729         v1ss=v1ss*whpb/wsc
1730         v2ss=v2ss*whpb/wsc
1731         v3ss=v3ss*whpb/wsc
1732       else
1733         ss_depth=ebr/whpb-0.25*eps(1,1)*wsc/whpb
1734       endif
1735
1736 C      if (lprint) then
1737       write (iout,'(/a)') "Disulfide bridge parameters:"
1738       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1739       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1740       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1741       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1742      & ' v3ss:',v3ss
1743 C      endif
1744       if (shield_mode.gt.0) then
1745       pi=3.141592d0
1746 C VSolvSphere the volume of solving sphere
1747 C      print *,pi,"pi"
1748 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
1749 C there will be no distinction between proline peptide group and normal peptide
1750 C group in case of shielding parameters
1751       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
1752       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
1753       write (iout,*) "VSolvSphere",VSolvSphere," VSolvSphere_div",
1754      &  VSolvSphere_div
1755 C long axis of side chain
1756       do i=1,ntyp
1757       long_r_sidechain(i)=vbldsc0(1,i)
1758       short_r_sidechain(i)=sigma0(i)
1759       enddo
1760       buff_shield=1.0d0
1761       endif
1762       return
1763   111 write (iout,*) "Error reading bending energy parameters."
1764       goto 999
1765   112 write (iout,*) "Error reading rotamer energy parameters."
1766       goto 999
1767   113 write (iout,*) "Error reading torsional energy parameters."
1768       goto 999
1769   114 write (iout,*) "Error reading double torsional energy parameters."
1770       goto 999
1771   115 write (iout,*)
1772      &  "Error reading cumulant (multibody energy) parameters."
1773       goto 999
1774   116 write (iout,*) "Error reading electrostatic energy parameters."
1775       goto 999
1776  1161 write (iout,*) "Error reading lipid parameters."
1777       goto 999
1778   117 write (iout,*) "Error reading side chain interaction parameters."
1779       goto 999
1780   118 write (iout,*) "Error reading SCp interaction parameters."
1781       goto 999
1782   119 write (iout,*) "Error reading SCCOR parameters"
1783       goto 999
1784   121 write (iout,*) "Error reading bond parameters"
1785   999 continue
1786 #ifdef MPI
1787       call MPI_Finalize(Ierror)
1788 #endif
1789       stop
1790       return
1791       end