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