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