source/wham/io_wham.F90

   1       module io_wham
   2
   3       use io_units
   4       use io_base
   5       use wham_data
   6 #ifndef CLUSTER
   7       use w_compar_data
   8 #endif
   9 !      use geometry_data
  10 !      use geometry
  11       implicit none
  12 !-----------------------------------------------------------------------------
  13 !
  14 !
  15 !-----------------------------------------------------------------------------
  16       contains
  17 !-----------------------------------------------------------------------------
  18 ! openunits.F
  19 !-----------------------------------------------------------------------------
  20 #ifndef CLUSTER
  21       subroutine openunits
  22 #ifdef WIN
  23       use dfport
  24 #endif
  25 !      implicit real*8 (a-h,o-z)
  26 !      include 'DIMENSIONS'
  27 !      include 'DIMENSIONS.ZSCOPT'
  28 #ifdef MPI
  29       use MPI_data
  30       include 'mpif.h'
  31 !      include 'COMMON.MPI'
  32 !      integer :: MyRank
  33       character(len=3) :: liczba
  34 #endif
  35 !      include 'COMMON.IOUNITS'
  36       integer :: lenpre,lenpot !,ilen
  37 !el      external ilen
  38
  39 #ifdef MPI
  40       MyRank=Me
  41 #endif
  42       call mygetenv('PREFIX',prefix)
  43       call mygetenv('SCRATCHDIR',scratchdir)
  44       call mygetenv('POT',pot)
  45       lenpre=ilen(prefix)
  46       lenpot=ilen(pot)
  47       call mygetenv('POT',pot)
  48       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
  49 ! Get the names and open the input files
  50       open (1,file=prefix(:ilen(prefix))//'.inp',status='old')
  51 ! Get parameter filenames and open the parameter files.
  52       call mygetenv('BONDPAR',bondname)
  53       open (ibond,file=bondname,status='old')
  54       call mygetenv('THETPAR',thetname)
  55       open (ithep,file=thetname,status='old')
  56       call mygetenv('ROTPAR',rotname)
  57       open (irotam,file=rotname,status='old')
  58       call mygetenv('TORPAR',torname)
  59       open (itorp,file=torname,status='old')
  60       call mygetenv('TORDPAR',tordname)
  61       open (itordp,file=tordname,status='old')
  62       call mygetenv('FOURIER',fouriername)
  63       open (ifourier,file=fouriername,status='old')
  64       call mygetenv('SCCORPAR',sccorname)
  65       open (isccor,file=sccorname,status='old')
  66       call mygetenv('ELEPAR',elename)
  67       open (ielep,file=elename,status='old')
  68       call mygetenv('SIDEPAR',sidename)
  69       open (isidep,file=sidename,status='old')
  70       call mygetenv('SIDEP',sidepname)
  71       open (isidep1,file=sidepname,status="old")
  72 #ifndef OLDSCP
  73 !
  74 ! 8/9/01 In the newest version SCp interaction constants are read from a file
  75 ! Use -DOLDSCP to use hard-coded constants instead.
  76 !
  77       call mygetenv('SCPPAR',scpname)
  78       open (iscpp,file=scpname,status='old')
  79 #endif
  80 #ifdef MPL
  81       if (MyID.eq.BossID) then
  82       MyRank = MyID/fgProcs
  83 #endif
  84 #ifdef MPI
  85       print *,'OpenUnits: processor',MyRank
  86       call numstr(MyRank,liczba)
  87       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//liczba
  88 #else
  89       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
  90 #endif
  91       open(iout,file=outname,status='unknown')
  92       write (iout,'(80(1h-))')
  93       write (iout,'(30x,a)') "FILE ASSIGNMENT"
  94       write (iout,'(80(1h-))')
  95       write (iout,*) "Input file                      : ",&
  96         prefix(:ilen(prefix))//'.inp'
  97       write (iout,*) "Output file                     : ",&
  98         outname(:ilen(outname))
  99       write (iout,*)
 100       write (iout,*) "Sidechain potential file        : ",&
 101         sidename(:ilen(sidename))
 102 #ifndef OLDSCP
 103       write (iout,*) "SCp potential file              : ",&
 104         scpname(:ilen(scpname))
 105 #endif
 106       write (iout,*) "Electrostatic potential file    : ",&
 107         elename(:ilen(elename))
 108       write (iout,*) "Cumulant coefficient file       : ",&
 109         fouriername(:ilen(fouriername))
 110       write (iout,*) "Torsional parameter file        : ",&
 111         torname(:ilen(torname))
 112       write (iout,*) "Double torsional parameter file : ",&
 113         tordname(:ilen(tordname))
 114       write (iout,*) "Backbone-rotamer parameter file : ",&
 115         sccorname(:ilen(sccorname))
 116       write (iout,*) "Bond & inertia constant file    : ",&
 117         bondname(:ilen(bondname))
 118       write (iout,*) "Bending parameter file          : ",&
 119         thetname(:ilen(thetname))
 120       write (iout,*) "Rotamer parameter file          : ",&
 121         rotname(:ilen(rotname))
 122       write (iout,'(80(1h-))')
 123       write (iout,*)
 124       return
 125       end subroutine openunits
 126 !-----------------------------------------------------------------------------
 127 ! molread_zs.F
 128 !-----------------------------------------------------------------------------
 129       subroutine molread(*)
 130 !
 131 ! Read molecular data.
 132 !
 133       use energy_data
 134       use geometry_data, only:nres,deg2rad,c,dc
 135       use control_data, only:iscode
 136       use control, only:rescode,setup_var,init_int_table
 137       use geometry, only:alloc_geo_arrays
 138       use energy, only:alloc_ener_arrays
 139 !      implicit real*8 (a-h,o-z)
 140 !      include 'DIMENSIONS'
 141 !      include 'DIMENSIONS.ZSCOPT'
 142 !      include 'COMMON.IOUNITS'
 143 !      include 'COMMON.GEO'
 144 !      include 'COMMON.VAR'
 145 !      include 'COMMON.INTERACT'
 146 !      include 'COMMON.LOCAL'
 147 !      include 'COMMON.NAMES'
 148 !      include 'COMMON.CHAIN'
 149 !      include 'COMMON.FFIELD'
 150 !      include 'COMMON.SBRIDGE'
 151 !      include 'COMMON.TORCNSTR'
 152 !      include 'COMMON.CONTROL'
 153       character(len=4),dimension(:),allocatable :: sequence !(nres)
 154 !el      integer :: rescode
 155 !el      real(kind=8) :: x(maxvar)
 156       character(len=320) :: controlcard !,ucase
 157       integer,dimension(nres) :: itype_pdb !(maxres)
 158       integer :: i,j,i1,i2,it1,it2
 159       real(kind=8) :: scalscp
 160 !el      logical :: seq_comp
 161       call card_concat(controlcard,.true.)
 162       call reada(controlcard,'SCAL14',scal14,0.4d0)
 163       call reada(controlcard,'SCALSCP',scalscp,1.0d0)
 164       call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
 165       call reada(controlcard,'TEMP0',temp0,300.0d0) !el
 166       call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
 167       r0_corr=cutoff_corr-delt_corr
 168       call readi(controlcard,"NRES",nres,0)
 169       allocate(sequence(nres+1))
 170 !el znamy juz ilosc reszt wiec mozna zaalokowac tablice do liczenia enerii
 171       call alloc_geo_arrays
 172       call alloc_ener_arrays
 173 ! alokacja dodatkowych tablic, ktore w unresie byly alokowanie w locie
 174 !----------------------------
 175       allocate(c(3,2*nres+2))
 176       allocate(dc(3,0:2*nres+2))
 177       allocate(itype(nres+2))
 178       allocate(itel(nres+2))
 179 !
 180 ! Zero out tableis.
 181       do i=1,2*nres+2
 182         do j=1,3
 183           c(j,i)=0.0D0
 184           dc(j,i)=0.0D0
 185         enddo
 186       enddo
 187       do i=1,nres+2
 188         itype(i)=0
 189         itel(i)=0
 190       enddo
 191 !--------------------------
 192 !
 193       iscode=index(controlcard,"ONE_LETTER")
 194       if (nres.le.0) then
 195         write (iout,*) "Error: no residues in molecule"
 196         return 1
 197       endif
 198       if (nres.gt.maxres) then
 199         write (iout,*) "Error: too many residues",nres,maxres
 200       endif
 201       write(iout,*) 'nres=',nres
 202 ! Read sequence of the protein
 203       if (iscode.gt.0) then
 204         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 205       else
 206         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 207       endif
 208 ! Convert sequence to numeric code
 209       do i=1,nres
 210         itype(i)=rescode(i,sequence(i),iscode)
 211       enddo
 212       write (iout,*) "Numeric code:"
 213       write (iout,'(20i4)') (itype(i),i=1,nres)
 214       do i=1,nres-1
 215 #ifdef PROCOR
 216         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 217 #else
 218         if (itype(i).eq.ntyp1) then
 219 #endif
 220           itel(i)=0
 221 #ifdef PROCOR
 222         else if (iabs(itype(i+1)).ne.20) then
 223 #else
 224         else if (iabs(itype(i)).ne.20) then
 225 #endif
 226           itel(i)=1
 227         else
 228           itel(i)=2
 229         endif
 230       enddo
 231        write (iout,*) "ITEL"
 232        do i=1,nres-1
 233          write (iout,*) i,itype(i),itel(i)
 234        enddo
 235       call read_bridge
 236
 237       if (with_dihed_constr) then
 238
 239       read (inp,*) ndih_constr
 240       if (ndih_constr.gt.0) then
 241         read (inp,*) ftors
 242         write (iout,*) 'FTORS',ftors
 243         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 244         write (iout,*) &
 245          'There are',ndih_constr,' constraints on phi angles.'
 246         do i=1,ndih_constr
 247           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 248         enddo
 249         do i=1,ndih_constr
 250           phi0(i)=deg2rad*phi0(i)
 251           drange(i)=deg2rad*drange(i)
 252         enddo
 253       endif
 254
 255       endif
 256
 257       nnt=1
 258       nct=nres
 259       if (itype(1).eq.ntyp1) nnt=2
 260       if (itype(nres).eq.ntyp1) nct=nct-1
 261       write(iout,*) 'NNT=',NNT,' NCT=',NCT
 262       call setup_var
 263       call init_int_table
 264       if (ns.gt.0) then
 265         write (iout,'(/a,i3,a)') 'The chain contains',ns,&
 266         ' disulfide-bridging cysteines.'
 267         write (iout,'(20i4)') (iss(i),i=1,ns)
 268         write (iout,'(/a/)') 'Pre-formed links are:'
 269         do i=1,nss
 270           i1=ihpb(i)-nres
 271           i2=jhpb(i)-nres
 272           it1=itype(i1)
 273           it2=itype(i2)
 274          write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
 275           restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',&
 276           dhpb(i),ebr,forcon(i)
 277         enddo
 278       endif
 279       write (iout,'(a)')
 280       return
 281       end subroutine molread
 282 !-----------------------------------------------------------------------------
 283 ! parmread.F
 284 !-----------------------------------------------------------------------------
 285       subroutine parmread(iparm,*)
 286 #else
 287       subroutine parmread
 288 #endif
 289 !
 290 ! Read the parameters of the probability distributions of the virtual-bond
 291 ! valence angles and the side chains and energy parameters.
 292 !
 293       use wham_data
 294
 295       use geometry_data
 296       use energy_data
 297       use control_data, only: maxtor,maxterm,maxlor,maxterm_sccor,&
 298           maxtermd_1,maxtermd_2,maxthetyp,maxthetyp1
 299       use MD_data
 300 !el      use MPI_data
 301 !el      use map_data
 302       use io_config, only: printmat
 303       use control, only: getenv_loc
 304
 305 #ifdef MPI
 306       use MPI_data
 307       include "mpif.h"
 308       integer :: IERROR
 309 #endif
 310 !      implicit real*8 (a-h,o-z)
 311 !      include 'DIMENSIONS'
 312 !      include 'DIMENSIONS.ZSCOPT'
 313 !      include 'DIMENSIONS.FREE'
 314 !      include 'COMMON.IOUNITS'
 315 !      include 'COMMON.CHAIN'
 316 !      include 'COMMON.INTERACT'
 317 !      include 'COMMON.GEO'
 318 !      include 'COMMON.LOCAL'
 319 !      include 'COMMON.TORSION'
 320 !      include 'COMMON.FFIELD'
 321 !      include 'COMMON.NAMES'
 322 !      include 'COMMON.SBRIDGE'
 323 !      include 'COMMON.WEIGHTS'
 324 !      include 'COMMON.ENEPS'
 325 !      include 'COMMON.SCCOR'
 326 !      include 'COMMON.SCROT'
 327 !      include 'COMMON.FREE'
 328       character(len=1) :: t1,t2,t3
 329       character(len=1) :: onelett(4) = (/"G","A","P","D"/)
 330       character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
 331       logical :: lprint
 332       real(kind=8),dimension(3,3,maxlob) :: blower      !(3,3,maxlob)
 333       character(len=800) :: controlcard
 334       character(len=256) :: bondname_t,thetname_t,rotname_t,torname_t,&
 335         tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,&
 336         sccorname_t
 337 !el      integer ilen
 338 !el   external ilen
 339       character(len=16) :: key
 340       integer :: iparm
 341 !el      real(kind=8) :: ip,mp
 342       real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,&
 343                 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm
 344       real(kind=8) :: v0ij,v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,rjunk,&
 345                 res1
 346       integer :: i,j,ichir1,ichir2,k,l,m,kk,ii,mm,junk,lll,ll,llll,n
 347       integer :: nlobi,iblock,maxinter,iscprol
 348 !
 349 ! Body
 350 !
 351 ! Set LPRINT=.TRUE. for debugging
 352       dwa16=2.0d0**(1.0d0/6.0d0)
 353       lprint=.false.
 354       itypro=20
 355 ! Assign virtual-bond length
 356       vbl=3.8D0
 357       vblinv=1.0D0/vbl
 358       vblinv2=vblinv*vblinv
 359 #ifndef CLUSTER
 360       call card_concat(controlcard,.true.)
 361       wname(4)="WCORRH"
 362 !el
 363 allocate(ww(max_eneW))
 364       do i=1,n_eneW
 365         key = wname(i)(:ilen(wname(i)))
 366         call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
 367       enddo
 368
 369       write (iout,*) "iparm",iparm," myparm",myparm
 370 ! If reading not own parameters, skip assignment
 371
 372       if (iparm.eq.myparm .or. .not.separate_parset) then
 373
 374 !
 375 ! Setup weights for UNRES
 376 !
 377       wsc=ww(1)
 378       wscp=ww(2)
 379       welec=ww(3)
 380       wcorr=ww(4)
 381       wcorr5=ww(5)
 382       wcorr6=ww(6)
 383       wel_loc=ww(7)
 384       wturn3=ww(8)
 385       wturn4=ww(9)
 386       wturn6=ww(10)
 387       wang=ww(11)
 388       wscloc=ww(12)
 389       wtor=ww(13)
 390       wtor_d=ww(14)
 391       wvdwpp=ww(16)
 392       wbond=ww(18)
 393       wsccor=ww(19)
 394
 395       endif
 396 !
 397 !el------
 398       allocate(weights(n_ene))
 399       weights(1)=wsc
 400       weights(2)=wscp
 401       weights(3)=welec
 402       weights(4)=wcorr
 403       weights(5)=wcorr5
 404       weights(6)=wcorr6
 405       weights(7)=wel_loc
 406       weights(8)=wturn3
 407       weights(9)=wturn4
 408       weights(10)=wturn6
 409       weights(11)=wang
 410       weights(12)=wscloc
 411       weights(13)=wtor
 412       weights(14)=wtor_d
 413       weights(15)=0 !wstrain !
 414       weights(16)=0 !wvdwpp !
 415       weights(17)=wbond
 416       weights(18)=0 !scal14 !
 417       weights(21)=wsccor
 418 ! el--------
 419       call card_concat(controlcard,.false.)
 420
 421 ! Return if not own parameters
 422
 423       if (iparm.ne.myparm .and. separate_parset) return
 424
 425       call reads(controlcard,"BONDPAR",bondname_t,bondname)
 426       open (ibond,file=bondname_t,status='old')
 427       rewind(ibond)
 428       call reads(controlcard,"THETPAR",thetname_t,thetname)
 429       open (ithep,file=thetname_t,status='old')
 430       rewind(ithep)
 431       call reads(controlcard,"ROTPAR",rotname_t,rotname)
 432       open (irotam,file=rotname_t,status='old')
 433       rewind(irotam)
 434       call reads(controlcard,"TORPAR",torname_t,torname)
 435       open (itorp,file=torname_t,status='old')
 436       rewind(itorp)
 437       call reads(controlcard,"TORDPAR",tordname_t,tordname)
 438       open (itordp,file=tordname_t,status='old')
 439       rewind(itordp)
 440       call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
 441       open (isccor,file=sccorname_t,status='old')
 442       rewind(isccor)
 443       call reads(controlcard,"FOURIER",fouriername_t,fouriername)
 444       open (ifourier,file=fouriername_t,status='old')
 445       rewind(ifourier)
 446       call reads(controlcard,"ELEPAR",elename_t,elename)
 447       open (ielep,file=elename_t,status='old')
 448       rewind(ielep)
 449       call reads(controlcard,"SIDEPAR",sidename_t,sidename)
 450       open (isidep,file=sidename_t,status='old')
 451       rewind(isidep)
 452       call reads(controlcard,"SCPPAR",scpname_t,scpname)
 453       open (iscpp,file=scpname_t,status='old')
 454       rewind(iscpp)
 455       write (iout,*) "Parameter set:",iparm
 456       write (iout,*) "Energy-term weights:"
 457       do i=1,n_eneW
 458         write (iout,'(a16,f10.5)') wname(i),ww(i)
 459       enddo
 460       write (iout,*) "Sidechain potential file        : ",&
 461         sidename_t(:ilen(sidename_t))
 462 #ifndef OLDSCP
 463       write (iout,*) "SCp potential file              : ",&
 464         scpname_t(:ilen(scpname_t))
 465 #endif
 466       write (iout,*) "Electrostatic potential file    : ",&
 467         elename_t(:ilen(elename_t))
 468       write (iout,*) "Cumulant coefficient file       : ",&
 469         fouriername_t(:ilen(fouriername_t))
 470       write (iout,*) "Torsional parameter file        : ",&
 471         torname_t(:ilen(torname_t))
 472       write (iout,*) "Double torsional parameter file : ",&
 473         tordname_t(:ilen(tordname_t))
 474       write (iout,*) "Backbone-rotamer parameter file : ",&
 475         sccorname(:ilen(sccorname))
 476       write (iout,*) "Bond & inertia constant file    : ",&
 477         bondname_t(:ilen(bondname_t))
 478       write (iout,*) "Bending parameter file          : ",&
 479         thetname_t(:ilen(thetname_t))
 480       write (iout,*) "Rotamer parameter file          : ",&
 481         rotname_t(:ilen(rotname_t))
 482 #endif
 483 !
 484 ! Read the virtual-bond parameters, masses, and moments of inertia
 485 ! and Stokes' radii of the peptide group and side chains
 486 !
 487       allocate(dsc(ntyp1)) !(ntyp1)
 488       allocate(dsc_inv(ntyp1)) !(ntyp1)
 489       allocate(nbondterm(ntyp)) !(ntyp)
 490       allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 491       allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 492 !el      allocate(msc(ntyp+1)) !(ntyp+1)
 493 !el      allocate(isc(ntyp+1)) !(ntyp+1)
 494 !el      allocate(restok(ntyp+1)) !(ntyp+1)
 495       allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 496
 497 #ifdef CRYST_BOND
 498       read (ibond,*) vbldp0,akp
 499       do i=1,ntyp
 500         nbondterm(i)=1
 501         read (ibond,*) vbldsc0(1,i),aksc(1,i)
 502         dsc(i) = vbldsc0(1,i)
 503         if (i.eq.10) then
 504           dsc_inv(i)=0.0D0
 505         else
 506           dsc_inv(i)=1.0D0/dsc(i)
 507         endif
 508       enddo
 509 #else
 510       read (ibond,*) ijunk,vbldp0,akp,rjunk
 511       do i=1,ntyp
 512         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
 513          j=1,nbondterm(i))
 514         dsc(i) = vbldsc0(1,i)
 515         if (i.eq.10) then
 516           dsc_inv(i)=0.0D0
 517         else
 518           dsc_inv(i)=1.0D0/dsc(i)
 519         endif
 520       enddo
 521 #endif
 522       if (lprint) then
 523         write(iout,'(/a/)')"Force constants virtual bonds:"
 524         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',&
 525          'inertia','Pstok'
 526         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
 527         do i=1,ntyp
 528           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),&
 529             vbldsc0(1,i),aksc(1,i),abond0(1,i)
 530           do j=2,nbondterm(i)
 531             write (iout,'(13x,3f10.5)') &
 532               vbldsc0(j,i),aksc(j,i),abond0(j,i)
 533           enddo
 534         enddo
 535       endif
 536 !----------------------------------------------------
 537       allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
 538       allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp))      !(-ntyp:ntyp)
 539       allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
 540       allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
 541       allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
 542       allocate(gthet(3,-ntyp:ntyp))     !(3,-ntyp:ntyp)
 543       do i=-ntyp,ntyp
 544         a0thet(i)=0.0D0
 545         do j=1,2
 546          do ichir1=-1,1
 547           do ichir2=-1,1
 548           athet(j,i,ichir1,ichir2)=0.0D0
 549           bthet(j,i,ichir1,ichir2)=0.0D0
 550           enddo
 551          enddo
 552         enddo
 553         do j=0,3
 554           polthet(j,i)=0.0D0
 555         enddo
 556         do j=1,3
 557           gthet(j,i)=0.0D0
 558         enddo
 559         theta0(i)=0.0D0
 560         sig0(i)=0.0D0
 561         sigc0(i)=0.0D0
 562       enddo
 563 !elwrite(iout,*) "parmread kontrol"
 564
 565 #ifdef CRYST_THETA
 566 !
 567 ! Read the parameters of the probability distribution/energy expression
 568 ! of the virtual-bond valence angles theta
 569 !
 570       do i=1,ntyp
 571         read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),&
 572           (bthet(j,i,1,1),j=1,2)
 573         read (ithep,*) (polthet(j,i),j=0,3)
 574 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 575         read (ithep,*) (gthet(j,i),j=1,3)
 576 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 577         read (ithep,*) theta0(i),sig0(i),sigc0(i)
 578         sigc0(i)=sigc0(i)**2
 579 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 580       enddo
 581 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 582       do i=1,ntyp
 583       athet(1,i,1,-1)=athet(1,i,1,1)
 584       athet(2,i,1,-1)=athet(2,i,1,1)
 585       bthet(1,i,1,-1)=-bthet(1,i,1,1)
 586       bthet(2,i,1,-1)=-bthet(2,i,1,1)
 587       athet(1,i,-1,1)=-athet(1,i,1,1)
 588       athet(2,i,-1,1)=-athet(2,i,1,1)
 589       bthet(1,i,-1,1)=bthet(1,i,1,1)
 590       bthet(2,i,-1,1)=bthet(2,i,1,1)
 591       enddo
 592 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 593       do i=-ntyp,-1
 594       a0thet(i)=a0thet(-i)
 595       athet(1,i,-1,-1)=athet(1,-i,1,1)
 596       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 597       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 598       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 599       athet(1,i,-1,1)=athet(1,-i,1,1)
 600       athet(2,i,-1,1)=-athet(2,-i,1,1)
 601       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 602       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 603       athet(1,i,1,-1)=-athet(1,-i,1,1)
 604       athet(2,i,1,-1)=athet(2,-i,1,1)
 605       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 606       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 607       theta0(i)=theta0(-i)
 608       sig0(i)=sig0(-i)
 609       sigc0(i)=sigc0(-i)
 610        do j=0,3
 611         polthet(j,i)=polthet(j,-i)
 612        enddo
 613        do j=1,3
 614          gthet(j,i)=gthet(j,-i)
 615        enddo
 616       enddo
 617 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 618       close (ithep)
 619 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 620       if (lprint) then
 621 !       write (iout,'(a)')
 622 !    &    'Parameters of the virtual-bond valence angles:'
 623 !       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
 624 !    & '    ATHETA0   ','         A1   ','        A2    ',
 625 !    & '        B1    ','         B2   '
 626 !       do i=1,ntyp
 627 !         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 628 !    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
 629 !       enddo
 630 !       write (iout,'(/a/9x,5a/79(1h-))')
 631 !    & 'Parameters of the expression for sigma(theta_c):',
 632 !    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
 633 !    & '     ALPH3    ','    SIGMA0C   '
 634 !       do i=1,ntyp
 635 !         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 636 !    &      (polthet(j,i),j=0,3),sigc0(i)
 637 !       enddo
 638 !       write (iout,'(/a/9x,5a/79(1h-))')
 639 !    & 'Parameters of the second gaussian:',
 640 !    & '    THETA0    ','     SIGMA0   ','        G1    ',
 641 !    & '        G2    ','         G3   '
 642 !       do i=1,ntyp
 643 !         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
 644 !    &       sig0(i),(gthet(j,i),j=1,3)
 645 !       enddo
 646         write (iout,'(a)') &
 647           'Parameters of the virtual-bond valence angles:'
 648         write (iout,'(/a/9x,5a/79(1h-))') &
 649        'Coefficients of expansion',&
 650        '     theta0   ','    a1*10^2   ','   a2*10^2    ',&
 651        '   b1*10^1    ','    b2*10^1   '
 652         do i=1,ntyp
 653           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
 654               a0thet(i),(100*athet(j,i,1,1),j=1,2),&
 655               (10*bthet(j,i,1,1),j=1,2)
 656         enddo
 657         write (iout,'(/a/9x,5a/79(1h-))') &
 658        'Parameters of the expression for sigma(theta_c):',&
 659        ' alpha0       ','  alph1       ',' alph2        ',&
 660        ' alhp3        ','   sigma0c    '
 661         do i=1,ntyp
 662           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
 663             (polthet(j,i),j=0,3),sigc0(i)
 664         enddo
 665         write (iout,'(/a/9x,5a/79(1h-))') &
 666        'Parameters of the second gaussian:',&
 667        '    theta0    ','  sigma0*10^2 ','      G1*10^-1',&
 668        '        G2    ','   G3*10^1    '
 669         do i=1,ntyp
 670           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
 671              100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 672         enddo
 673       endif
 674 #else
 675 !
 676 ! Read the parameters of Utheta determined from ab initio surfaces
 677 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 678 !
 679 !      write (iout,*) "tu dochodze"
 680       read (ithep,*) nthetyp,ntheterm,ntheterm2,&
 681         ntheterm3,nsingle,ndouble
 682       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 683
 684 !----------------------------------------------------
 685       allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
 686       allocate(aa0thet(-maxthetyp1:maxthetyp1,&
 687         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 688 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 689       allocate(aathet(ntheterm,-maxthetyp1:maxthetyp1,&
 690         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 691 !(maxtheterm,-maxthetyp1:maxthetyp1,&
 692 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 693       allocate(bbthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
 694         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 695       allocate(ccthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
 696         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 697       allocate(ddthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
 698         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 699       allocate(eethet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
 700         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 701 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
 702 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 703       allocate(ffthet(ndouble,ndouble,ntheterm3,-maxthetyp1:maxthetyp1,&
 704         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 705       allocate(ggthet(ndouble,ndouble,ntheterm3,-maxthetyp1:maxthetyp1,&
 706         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 707 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
 708 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 709
 710
 711       read (ithep,*) (ithetyp(i),i=1,ntyp1)
 712       do i=-ntyp1,-1
 713         ithetyp(i)=-ithetyp(-i)
 714       enddo
 715 !      write (iout,*) "tu dochodze"
 716       do iblock=1,2
 717       do i=-maxthetyp,maxthetyp
 718         do j=-maxthetyp,maxthetyp
 719           do k=-maxthetyp,maxthetyp
 720             aa0thet(i,j,k,iblock)=0.0d0
 721             do l=1,ntheterm
 722               aathet(l,i,j,k,iblock)=0.0d0
 723             enddo
 724             do l=1,ntheterm2
 725               do m=1,nsingle
 726                 bbthet(m,l,i,j,k,iblock)=0.0d0
 727                 ccthet(m,l,i,j,k,iblock)=0.0d0
 728                 ddthet(m,l,i,j,k,iblock)=0.0d0
 729                 eethet(m,l,i,j,k,iblock)=0.0d0
 730               enddo
 731             enddo
 732             do l=1,ntheterm3
 733               do m=1,ndouble
 734                 do mm=1,ndouble
 735                  ffthet(mm,m,l,i,j,k,iblock)=0.0d0
 736                  ggthet(mm,m,l,i,j,k,iblock)=0.0d0
 737                 enddo
 738               enddo
 739             enddo
 740           enddo
 741         enddo
 742       enddo
 743       enddo
 744       do iblock=1,2
 745       do i=0,nthetyp
 746         do j=-nthetyp,nthetyp
 747           do k=-nthetyp,nthetyp
 748             read (ithep,'(6a)') res1
 749             read (ithep,*) aa0thet(i,j,k,iblock)
 750             read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
 751             read (ithep,*) &
 752              ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
 753               (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
 754               (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
 755               (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
 756               ll=1,ntheterm2)
 757             read (ithep,*) &
 758             (((ffthet(llll,lll,ll,i,j,k,iblock),&
 759                ffthet(lll,llll,ll,i,j,k,iblock),&
 760                ggthet(llll,lll,ll,i,j,k,iblock),&
 761                ggthet(lll,llll,ll,i,j,k,iblock),&
 762                llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 763           enddo
 764         enddo
 765       enddo
 766 !
 767 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
 768 ! coefficients of theta-and-gamma-dependent terms are zero.
 769 !
 770       do i=1,nthetyp
 771         do j=1,nthetyp
 772           do l=1,ntheterm
 773             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
 774             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
 775           enddo
 776           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
 777           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
 778         enddo
 779         do l=1,ntheterm
 780           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 781         enddo
 782         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 783       enddo
 784       enddo
 785 ! Substitution for D aminoacids from symmetry.
 786       do iblock=1,2
 787       do i=-nthetyp,0
 788         do j=-nthetyp,nthetyp
 789           do k=-nthetyp,nthetyp
 790            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
 791            do l=1,ntheterm
 792            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
 793            enddo
 794            do ll=1,ntheterm2
 795             do lll=1,nsingle
 796             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
 797             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
 798             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
 799             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
 800             enddo
 801           enddo
 802           do ll=1,ntheterm3
 803            do lll=2,ndouble
 804             do llll=1,lll-1
 805             ffthet(llll,lll,ll,i,j,k,iblock)= &
 806             ffthet(llll,lll,ll,-i,-j,-k,iblock)
 807             ffthet(lll,llll,ll,i,j,k,iblock)= &
 808             ffthet(lll,llll,ll,-i,-j,-k,iblock)
 809             ggthet(llll,lll,ll,i,j,k,iblock)= &
 810             -ggthet(llll,lll,ll,-i,-j,-k,iblock)
 811             ggthet(lll,llll,ll,i,j,k,iblock)= &
 812             -ggthet(lll,llll,ll,-i,-j,-k,iblock)
 813             enddo !ll
 814            enddo  !lll
 815           enddo   !llll
 816          enddo    !k
 817         enddo     !j
 818        enddo      !i
 819       enddo       !iblock
 820
 821 !
 822 ! Control printout of the coefficients of virtual-bond-angle potentials
 823 !
 824 do iblock=1,2
 825       if (lprint) then
 826         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 827         do i=1,nthetyp+1
 828           do j=1,nthetyp+1
 829             do k=1,nthetyp+1
 830               write (iout,'(//4a)') &
 831                'Type ',onelett(i),onelett(j),onelett(k)
 832               write (iout,'(//a,10x,a)') " l","a[l]"
 833               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
 834               write (iout,'(i2,1pe15.5)') &
 835                  (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
 836             do l=1,ntheterm2
 837               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))') &
 838                 "b",l,"c",l,"d",l,"e",l
 839               do m=1,nsingle
 840                 write (iout,'(i2,4(1pe15.5))') m,&
 841                 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
 842                 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
 843               enddo
 844             enddo
 845             do l=1,ntheterm3
 846               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))') &
 847                 "f+",l,"f-",l,"g+",l,"g-",l
 848               do m=2,ndouble
 849                 do n=1,m-1
 850                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
 851                     ffthet(n,m,l,i,j,k,iblock),&
 852                     ffthet(m,n,l,i,j,k,iblock),&
 853                     ggthet(n,m,l,i,j,k,iblock),&
 854                     ggthet(m,n,l,i,j,k,iblock)
 855                 enddo
 856               enddo
 857             enddo
 858           enddo
 859         enddo
 860       enddo
 861       call flush(iout)
 862       endif
 863 enddo
 864 #endif
 865 !-------------------------------------------
 866       allocate(nlob(ntyp1)) !(ntyp1)
 867       allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
 868       allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
 869       allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
 870
 871       do i=1,ntyp
 872         do j=1,maxlob
 873           bsc(j,i)=0.0D0
 874           nlob(i)=0
 875         enddo
 876       enddo
 877       nlob(ntyp1)=0
 878       dsc(ntyp1)=0.0D0
 879
 880       do i=-ntyp,ntyp
 881         do j=1,maxlob
 882           do k=1,3
 883             censc(k,j,i)=0.0D0
 884           enddo
 885           do k=1,3
 886             do l=1,3
 887               gaussc(l,k,j,i)=0.0D0
 888             enddo
 889           enddo
 890         enddo
 891       enddo
 892
 893 #ifdef CRYST_SC
 894 !
 895 ! Read the parameters of the probability distribution/energy expression
 896 ! of the side chains.
 897 !
 898       do i=1,ntyp
 899 !c      write (iout,*) "tu dochodze",i
 900         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
 901         if (i.eq.10) then
 902           dsc_inv(i)=0.0D0
 903         else
 904           dsc_inv(i)=1.0D0/dsc(i)
 905         endif
 906         if (i.ne.10) then
 907         do j=1,nlob(i)
 908           do k=1,3
 909             do l=1,3
 910               blower(l,k,j)=0.0D0
 911             enddo
 912           enddo
 913         enddo
 914         bsc(1,i)=0.0D0
 915         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
 916         censc(1,1,-i)=censc(1,1,i)
 917         censc(2,1,-i)=censc(2,1,i)
 918         censc(3,1,-i)=-censc(3,1,i)
 919         do j=2,nlob(i)
 920           read (irotam,*) bsc(j,i)
 921           read (irotam,*) (censc(k,j,i),k=1,3),&
 922                                        ((blower(k,l,j),l=1,k),k=1,3)
 923         censc(1,j,-i)=censc(1,j,i)
 924         censc(2,j,-i)=censc(2,j,i)
 925         censc(3,j,-i)=-censc(3,j,i)
 926 ! BSC is amplitude of Gaussian
 927         enddo
 928         do j=1,nlob(i)
 929           do k=1,3
 930             do l=1,k
 931               akl=0.0D0
 932               do m=1,3
 933                 akl=akl+blower(k,m,j)*blower(l,m,j)
 934               enddo
 935               gaussc(k,l,j,i)=akl
 936               gaussc(l,k,j,i)=akl
 937              if (((k.eq.3).and.(l.ne.3)) &
 938               .or.((l.eq.3).and.(k.ne.3))) then
 939                 gaussc(k,l,j,-i)=-akl
 940                 gaussc(l,k,j,-i)=-akl
 941               else
 942                 gaussc(k,l,j,-i)=akl
 943                 gaussc(l,k,j,-i)=akl
 944               endif
 945             enddo
 946           enddo
 947         enddo
 948         endif
 949       enddo
 950       close (irotam)
 951       if (lprint) then
 952         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 953         do i=1,ntyp
 954           nlobi=nlob(i)
 955           if (nlobi.gt.0) then
 956           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
 957            ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 958 !          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 959 !          write (iout,'(a,f10.4,4(16x,f10.4))')
 960 !     &                             'Center  ',(bsc(j,i),j=1,nlobi)
 961 !          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
 962            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
 963                                    'log h',(bsc(j,i),j=1,nlobi)
 964            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
 965           'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 966 !          write (iout,'(a)')
 967 !         do j=1,nlobi
 968 !           ind=0
 969 !           do k=1,3
 970 !             do l=1,k
 971 !              ind=ind+1
 972 !              blower(k,l,j)=gaussc(ind,j,i)
 973 !             enddo
 974 !           enddo
 975 !         enddo
 976           do k=1,3
 977             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
 978                        ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 979           enddo
 980           endif
 981         enddo
 982       endif
 983 #else
 984 !
 985 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
 986 ! added by Urszula Kozlowska 07/11/2007
 987 !
 988       allocate(sc_parmin(65,ntyp))      !(maxsccoef,ntyp)
 989
 990       do i=1,ntyp
 991         read (irotam,*)
 992        if (i.eq.10) then
 993          read (irotam,*)
 994        else
 995          do j=1,65
 996            read(irotam,*) sc_parmin(j,i)
 997          enddo
 998        endif
 999       enddo
1000 #endif
1001       close(irotam)
1002 #ifdef CRYST_TOR
1003 !
1004 ! Read torsional parameters in old format
1005 !
1006       allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1007
1008       read (itorp,*) ntortyp,nterm_old
1009       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1010       read (itorp,*) (itortyp(i),i=1,ntyp)
1011
1012 !el from energy module--------
1013       allocate(v1(nterm_old,ntortyp,ntortyp))
1014       allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1015 !el---------------------------
1016
1017       do i=1,ntortyp
1018         do j=1,ntortyp
1019           read (itorp,'(a)')
1020           do k=1,nterm_old
1021             read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
1022           enddo
1023         enddo
1024       enddo
1025       close (itorp)
1026       if (lprint) then
1027         write (iout,'(/a/)') 'Torsional constants:'
1028         do i=1,ntortyp
1029           do j=1,ntortyp
1030             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1031             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1032           enddo
1033         enddo
1034       endif
1035
1036
1037 #else
1038 !
1039 ! Read torsional parameters
1040 !
1041       allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1042
1043       read (itorp,*) ntortyp
1044       read (itorp,*) (itortyp(i),i=1,ntyp)
1045       write (iout,*) 'ntortyp',ntortyp
1046
1047 !el from energy module---------
1048       allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1049       allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1050
1051       allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1052       allocate(vlor2(maxlor,ntortyp,ntortyp))
1053       allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1054       allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1055
1056       allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1057       allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1058 !el---------------------------
1059       do iblock=1,2
1060         do i=-ntortyp,ntortyp
1061           do j=-ntortyp,ntortyp
1062             nterm(i,j,iblock)=0
1063             nlor(i,j,iblock)=0
1064           enddo
1065         enddo
1066       enddo
1067 !el---------------------------
1068
1069       do iblock=1,2
1070       do i=-ntyp,-1
1071        itortyp(i)=-itortyp(-i)
1072       enddo
1073 !      write (iout,*) 'ntortyp',ntortyp
1074       do i=0,ntortyp-1
1075         do j=-ntortyp+1,ntortyp-1
1076           read (itorp,*) nterm(i,j,iblock),&
1077                 nlor(i,j,iblock)
1078           nterm(-i,-j,iblock)=nterm(i,j,iblock)
1079           nlor(-i,-j,iblock)=nlor(i,j,iblock)
1080           v0ij=0.0d0
1081           si=-1.0d0
1082           do k=1,nterm(i,j,iblock)
1083             read (itorp,*) kk,v1(k,i,j,iblock),&
1084             v2(k,i,j,iblock)
1085             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1086             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1087             v0ij=v0ij+si*v1(k,i,j,iblock)
1088             si=-si
1089          enddo
1090           do k=1,nlor(i,j,iblock)
1091             read (itorp,*) kk,vlor1(k,i,j),&
1092               vlor2(k,i,j),vlor3(k,i,j)
1093             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1094           enddo
1095           v0(i,j,iblock)=v0ij
1096           v0(-i,-j,iblock)=v0ij
1097         enddo
1098       enddo
1099       enddo
1100       close (itorp)
1101       if (lprint) then
1102         do iblock=1,2 !el
1103         write (iout,'(/a/)') 'Torsional constants:'
1104         do i=1,ntortyp
1105           do j=1,ntortyp
1106             write (iout,*) 'ityp',i,' jtyp',j
1107             write (iout,*) 'Fourier constants'
1108             do k=1,nterm(i,j,iblock)
1109               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1110               v2(k,i,j,iblock)
1111             enddo
1112             write (iout,*) 'Lorenz constants'
1113             do k=1,nlor(i,j,iblock)
1114               write (iout,'(3(1pe15.5))') &
1115                vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1116             enddo
1117           enddo
1118         enddo
1119         enddo
1120       endif
1121 !
1122 ! 6/23/01 Read parameters for double torsionals
1123 !
1124 !el from energy module------------
1125       allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1126       allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1127 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1128       allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1129       allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1130         !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1131       allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1132       allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1133         !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1134 !---------------------------------
1135
1136       do iblock=1,2
1137       do i=0,ntortyp-1
1138         do j=-ntortyp+1,ntortyp-1
1139           do k=-ntortyp+1,ntortyp-1
1140             read (itordp,'(3a1)') t1,t2,t3
1141 !              write (iout,*) "OK onelett",
1142 !     &         i,j,k,t1,t2,t3
1143
1144             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1145               .or. t3.ne.toronelet(k)) then
1146               write (iout,*) "Error in double torsional parameter file",&
1147                i,j,k,t1,t2,t3
1148 #ifdef MPI
1149               call MPI_Finalize(Ierror)
1150 #endif
1151                stop "Error in double torsional parameter file"
1152             endif
1153           read (itordp,*) ntermd_1(i,j,k,iblock),&
1154                ntermd_2(i,j,k,iblock)
1155             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1156             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1157             read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,&
1158                ntermd_1(i,j,k,iblock))
1159             read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,&
1160                ntermd_1(i,j,k,iblock))
1161             read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,&
1162                ntermd_1(i,j,k,iblock))
1163             read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,&
1164                ntermd_1(i,j,k,iblock))
1165 ! Martix of D parameters for one dimesional foureir series
1166             do l=1,ntermd_1(i,j,k,iblock)
1167              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1168              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1169              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1170              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1171 !            write(iout,*) "whcodze" ,
1172 !     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1173             enddo
1174             read (itordp,*) ((v2c(l,m,i,j,k,iblock),&
1175                v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1176                v2s(m,l,i,j,k,iblock),&
1177                m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1178 ! Martix of D parameters for two dimesional fourier series
1179             do l=1,ntermd_2(i,j,k,iblock)
1180              do m=1,l-1
1181              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1182              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1183              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1184              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1185              enddo!m
1186             enddo!l
1187           enddo!k
1188         enddo!j
1189       enddo!i
1190       enddo!iblock
1191       if (lprint) then
1192       write (iout,*)
1193       write (iout,*) 'Constants for double torsionals'
1194       do iblock=1,2
1195       do i=0,ntortyp-1
1196         do j=-ntortyp+1,ntortyp-1
1197           do k=-ntortyp+1,ntortyp-1
1198             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1199               ' nsingle',ntermd_1(i,j,k,iblock),&
1200               ' ndouble',ntermd_2(i,j,k,iblock)
1201             write (iout,*)
1202             write (iout,*) 'Single angles:'
1203             do l=1,ntermd_1(i,j,k,iblock)
1204               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1205                  v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1206                  v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1207                  v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1208             enddo
1209             write (iout,*)
1210             write (iout,*) 'Pairs of angles:'
1211             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1212             do l=1,ntermd_2(i,j,k,iblock)
1213               write (iout,'(i5,20f10.5)') &
1214                l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1215             enddo
1216             write (iout,*)
1217            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1218             do l=1,ntermd_2(i,j,k,iblock)
1219               write (iout,'(i5,20f10.5)') &
1220                l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1221                (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1222             enddo
1223             write (iout,*)
1224           enddo
1225         enddo
1226       enddo
1227       enddo
1228       endif
1229 #endif
1230 !elwrite(iout,*) "parmread kontrol sc-bb"
1231 ! Read of Side-chain backbone correlation parameters
1232 ! Modified 11 May 2012 by Adasko
1233 !CC
1234 !
1235      read (isccor,*) nsccortyp
1236
1237      maxinter=3
1238 !c maxinter is maximum interaction sites
1239 !write(iout,*)"maxterm_sccor",maxterm_sccor
1240 !el from module energy-------------
1241       allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1242       allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1243       allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1244       allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp))   !(maxterm_sccor,20,20)
1245 !-----------------------------------
1246       allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1247 !-----------------------------------
1248       allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1249       allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1250                -nsccortyp:nsccortyp))
1251       allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1252                -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1253       allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1254                -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1255 !-----------------------------------
1256       do i=-nsccortyp,nsccortyp
1257         do j=-nsccortyp,nsccortyp
1258           nterm_sccor(j,i)=0
1259         enddo
1260       enddo
1261 !-----------------------------------
1262
1263       read (isccor,*) (isccortyp(i),i=1,ntyp)
1264       do i=-ntyp,-1
1265         isccortyp(i)=-isccortyp(-i)
1266       enddo
1267       iscprol=isccortyp(20)
1268 !      write (iout,*) 'ntortyp',ntortyp
1269 !      maxinter=3
1270 !c maxinter is maximum interaction sites
1271       do l=1,maxinter
1272       do i=1,nsccortyp
1273         do j=1,nsccortyp
1274           read (isccor,*) &
1275       nterm_sccor(i,j),nlor_sccor(i,j)
1276           v0ijsccor=0.0d0
1277           v0ijsccor1=0.0d0
1278           v0ijsccor2=0.0d0
1279           v0ijsccor3=0.0d0
1280           si=-1.0d0
1281           nterm_sccor(-i,j)=nterm_sccor(i,j)
1282           nterm_sccor(-i,-j)=nterm_sccor(i,j)
1283           nterm_sccor(i,-j)=nterm_sccor(i,j)
1284           do k=1,nterm_sccor(i,j)
1285             read (isccor,*) kk,v1sccor(k,l,i,j),&
1286             v2sccor(k,l,i,j)
1287             if (j.eq.iscprol) then
1288              if (i.eq.isccortyp(10)) then
1289              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1290              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1291              else
1292              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1293                               +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1294              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1295                               +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1296              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1297              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1298              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1299              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1300              endif
1301             else
1302              if (i.eq.isccortyp(10)) then
1303              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1304              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1305              else
1306                if (j.eq.isccortyp(10)) then
1307              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1308              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1309                else
1310              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1311              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1312              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1313              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1314              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1315              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1316                 endif
1317                endif
1318             endif
1319             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1320             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1321             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1322             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1323             si=-si
1324            enddo
1325           do k=1,nlor_sccor(i,j)
1326             read (isccor,*) kk,vlor1sccor(k,i,j),&
1327               vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1328             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1329       (1+vlor3sccor(k,i,j)**2)
1330           enddo
1331           v0sccor(l,i,j)=v0ijsccor
1332           v0sccor(l,-i,j)=v0ijsccor1
1333           v0sccor(l,i,-j)=v0ijsccor2
1334           v0sccor(l,-i,-j)=v0ijsccor3
1335           enddo
1336         enddo
1337       enddo
1338       close (isccor)
1339       if (lprint) then
1340         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1341         do i=1,nsccortyp
1342           do j=1,nsccortyp
1343             write (iout,*) 'ityp',i,' jtyp',j
1344             write (iout,*) 'Fourier constants'
1345             do k=1,nterm_sccor(i,j)
1346               write (iout,'(2(1pe15.5))') &
1347          (v1sccor(k,l,i,j),v2sccor(k,l,i,j),l=1,maxinter)
1348             enddo
1349             write (iout,*) 'Lorenz constants'
1350             do k=1,nlor_sccor(i,j)
1351               write (iout,'(3(1pe15.5))') &
1352                vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1353             enddo
1354           enddo
1355         enddo
1356       endif
1357 !
1358 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1359 !         interaction energy of the Gly, Ala, and Pro prototypes.
1360 !
1361       read (ifourier,*) nloctyp
1362 !el write(iout,*)"nloctyp",nloctyp
1363 !el from module energy-------
1364       allocate(b1(2,-nloctyp-1:nloctyp+1))      !(2,-maxtor:maxtor)
1365       allocate(b2(2,-nloctyp-1:nloctyp+1))      !(2,-maxtor:maxtor)
1366       allocate(b1tilde(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1367       allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1368       allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1369       allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1370       allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1371       allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1372       do i=1,2
1373         do ii=-nloctyp-1,nloctyp+1
1374           b1(i,ii)=0.0d0
1375           b2(i,ii)=0.0d0
1376           b1tilde(i,ii)=0.0d0
1377           do j=1,2
1378             cc(j,i,ii)=0.0d0
1379             dd(j,i,ii)=0.0d0
1380             ee(j,i,ii)=0.0d0
1381             ctilde(j,i,ii)=0.0d0
1382             dtilde(j,i,ii)=0.0d0
1383           enddo
1384         enddo
1385       enddo
1386 !--------------------------------
1387       allocate(b(13,0:nloctyp))
1388
1389       do i=0,nloctyp-1
1390         read (ifourier,*)
1391         read (ifourier,*) (b(ii,i),ii=1,13)
1392         if (lprint) then
1393         write (iout,*) 'Type',i
1394         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1395         endif
1396         B1(1,i)  = b(3,i)
1397         B1(2,i)  = b(5,i)
1398         B1(1,-i) = b(3,i)
1399         B1(2,-i) = -b(5,i)
1400 !        b1(1,i)=0.0d0
1401 !        b1(2,i)=0.0d0
1402         B1tilde(1,i) = b(3,i)
1403         B1tilde(2,i) =-b(5,i)
1404         B1tilde(1,-i) =-b(3,i)
1405         B1tilde(2,-i) =b(5,i)
1406 !        b1tilde(1,i)=0.0d0
1407 !        b1tilde(2,i)=0.0d0
1408         B2(1,i)  = b(2,i)
1409         B2(2,i)  = b(4,i)
1410         B2(1,-i)  =b(2,i)
1411         B2(2,-i)  =-b(4,i)
1412
1413 !        b2(1,i)=0.0d0
1414 !        b2(2,i)=0.0d0
1415         CC(1,1,i)= b(7,i)
1416         CC(2,2,i)=-b(7,i)
1417         CC(2,1,i)= b(9,i)
1418         CC(1,2,i)= b(9,i)
1419         CC(1,1,-i)= b(7,i)
1420         CC(2,2,-i)=-b(7,i)
1421         CC(2,1,-i)=-b(9,i)
1422         CC(1,2,-i)=-b(9,i)
1423 !        CC(1,1,i)=0.0d0
1424 !        CC(2,2,i)=0.0d0
1425 !        CC(2,1,i)=0.0d0
1426 !        CC(1,2,i)=0.0d0
1427         Ctilde(1,1,i)=b(7,i)
1428         Ctilde(1,2,i)=b(9,i)
1429         Ctilde(2,1,i)=-b(9,i)
1430         Ctilde(2,2,i)=b(7,i)
1431         Ctilde(1,1,-i)=b(7,i)
1432         Ctilde(1,2,-i)=-b(9,i)
1433         Ctilde(2,1,-i)=b(9,i)
1434         Ctilde(2,2,-i)=b(7,i)
1435
1436 !        Ctilde(1,1,i)=0.0d0
1437 !        Ctilde(1,2,i)=0.0d0
1438 !        Ctilde(2,1,i)=0.0d0
1439 !        Ctilde(2,2,i)=0.0d0
1440         DD(1,1,i)= b(6,i)
1441         DD(2,2,i)=-b(6,i)
1442         DD(2,1,i)= b(8,i)
1443         DD(1,2,i)= b(8,i)
1444         DD(1,1,-i)= b(6,i)
1445         DD(2,2,-i)=-b(6,i)
1446         DD(2,1,-i)=-b(8,i)
1447         DD(1,2,-i)=-b(8,i)
1448 !        DD(1,1,i)=0.0d0
1449 !        DD(2,2,i)=0.0d0
1450 !        DD(2,1,i)=0.0d0
1451 !        DD(1,2,i)=0.0d0
1452         Dtilde(1,1,i)=b(6,i)
1453         Dtilde(1,2,i)=b(8,i)
1454         Dtilde(2,1,i)=-b(8,i)
1455         Dtilde(2,2,i)=b(6,i)
1456         Dtilde(1,1,-i)=b(6,i)
1457         Dtilde(1,2,-i)=-b(8,i)
1458         Dtilde(2,1,-i)=b(8,i)
1459         Dtilde(2,2,-i)=b(6,i)
1460
1461 !        Dtilde(1,1,i)=0.0d0
1462 !        Dtilde(1,2,i)=0.0d0
1463 !        Dtilde(2,1,i)=0.0d0
1464 !        Dtilde(2,2,i)=0.0d0
1465         EE(1,1,i)= b(10,i)+b(11,i)
1466         EE(2,2,i)=-b(10,i)+b(11,i)
1467         EE(2,1,i)= b(12,i)-b(13,i)
1468         EE(1,2,i)= b(12,i)+b(13,i)
1469         EE(1,1,-i)= b(10,i)+b(11,i)
1470         EE(2,2,-i)=-b(10,i)+b(11,i)
1471         EE(2,1,-i)=-b(12,i)+b(13,i)
1472         EE(1,2,-i)=-b(12,i)-b(13,i)
1473
1474 !        ee(1,1,i)=1.0d0
1475 !        ee(2,2,i)=1.0d0
1476 !        ee(2,1,i)=0.0d0
1477 !        ee(1,2,i)=0.0d0
1478 !        ee(2,1,i)=ee(1,2,i)
1479
1480       enddo
1481       if (lprint) then
1482       do i=1,nloctyp
1483         write (iout,*) 'Type',i
1484         write (iout,*) 'B1'
1485 !        write (iout,'(f10.5)') B1(:,i)
1486         write(iout,*) B1(1,i),B1(2,i)
1487         write (iout,*) 'B2'
1488 !        write (iout,'(f10.5)') B2(:,i)
1489         write(iout,*) B2(1,i),B2(2,i)
1490         write (iout,*) 'CC'
1491         do j=1,2
1492           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1493         enddo
1494         write(iout,*) 'DD'
1495         do j=1,2
1496           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1497         enddo
1498         write(iout,*) 'EE'
1499         do j=1,2
1500           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1501         enddo
1502       enddo
1503       endif
1504 !
1505 ! Read electrostatic-interaction parameters
1506 !
1507       if (lprint) then
1508         write (iout,'(/a)') 'Electrostatic interaction constants:'
1509         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
1510                   'IT','JT','APP','BPP','AEL6','AEL3'
1511       endif
1512       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
1513       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
1514       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
1515       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
1516       close (ielep)
1517       do i=1,2
1518         do j=1,2
1519         rri=rpp(i,j)**6
1520         app (i,j)=epp(i,j)*rri*rri
1521         bpp (i,j)=-2.0D0*epp(i,j)*rri
1522         ael6(i,j)=elpp6(i,j)*4.2D0**6
1523         ael3(i,j)=elpp3(i,j)*4.2D0**3
1524         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
1525                           ael6(i,j),ael3(i,j)
1526         enddo
1527       enddo
1528 !
1529 ! Read side-chain interaction parameters.
1530 !
1531 !el from module energy - COMMON.INTERACT-------
1532       allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
1533       allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
1534       allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
1535       allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
1536       allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
1537       do i=1,ntyp
1538         do j=1,ntyp
1539           augm(i,j)=0.0D0
1540         enddo
1541         chip(i)=0.0D0
1542         alp(i)=0.0D0
1543         sigma0(i)=0.0D0
1544         sigii(i)=0.0D0
1545         rr0(i)=0.0D0
1546       enddo
1547 !--------------------------------
1548
1549       read (isidep,*) ipot,expon
1550 !el      if (ipot.lt.1 .or. ipot.gt.5) then
1551 !        write (iout,'(2a)') 'Error while reading SC interaction',&
1552 !                     'potential file - unknown potential type.'
1553 !        stop
1554 !wl      endif
1555       expon2=expon/2
1556       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
1557        ', exponents are ',expon,2*expon
1558 !      goto (10,20,30,30,40) ipot
1559       select case(ipot)
1560 !----------------------- LJ potential ---------------------------------
1561        case (1)
1562 !   10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1563         read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1564         if (lprint) then
1565           write (iout,'(/a/)') 'Parameters of the LJ potential:'
1566           write (iout,'(a/)') 'The epsilon array:'
1567           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1568           write (iout,'(/a)') 'One-body parameters:'
1569           write (iout,'(a,4x,a)') 'residue','sigma'
1570           write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1571         endif
1572 !      goto 50
1573 !----------------------- LJK potential --------------------------------
1574        case (2)
1575 !   20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1576         read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1577           (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1578         if (lprint) then
1579           write (iout,'(/a/)') 'Parameters of the LJK potential:'
1580           write (iout,'(a/)') 'The epsilon array:'
1581           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1582           write (iout,'(/a)') 'One-body parameters:'
1583           write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1584           write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
1585                 i=1,ntyp)
1586         endif
1587 !      goto 50
1588 !---------------------- GB or BP potential -----------------------------
1589        case (3:4)
1590 !   30 do i=1,ntyp
1591         do i=1,ntyp
1592          read (isidep,*)(eps(i,j),j=i,ntyp)
1593         enddo
1594         read (isidep,*)(sigma0(i),i=1,ntyp)
1595         read (isidep,*)(sigii(i),i=1,ntyp)
1596         read (isidep,*)(chip(i),i=1,ntyp)
1597         read (isidep,*)(alp(i),i=1,ntyp)
1598 ! For the GB potential convert sigma'**2 into chi'
1599         if (ipot.eq.4) then
1600           do i=1,ntyp
1601             chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1602           enddo
1603         endif
1604         if (lprint) then
1605           write (iout,'(/a/)') 'Parameters of the BP potential:'
1606           write (iout,'(a/)') 'The epsilon array:'
1607           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1608           write (iout,'(/a)') 'One-body parameters:'
1609           write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',&
1610                '    chip  ','    alph  '
1611           write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
1612                            chip(i),alp(i),i=1,ntyp)
1613         endif
1614 !      goto 50
1615 !--------------------- GBV potential -----------------------------------
1616        case (5)
1617 !   40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1618         read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1619           (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
1620         (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1621         if (lprint) then
1622           write (iout,'(/a/)') 'Parameters of the GBV potential:'
1623           write (iout,'(a/)') 'The epsilon array:'
1624           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1625           write (iout,'(/a)') 'One-body parameters:'
1626           write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',&
1627             's||/s_|_^2','    chip  ','    alph  '
1628           write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
1629                  sigii(i),chip(i),alp(i),i=1,ntyp)
1630         endif
1631        case default
1632         write (iout,'(2a)') 'Error while reading SC interaction',&
1633                      'potential file - unknown potential type.'
1634         stop
1635 !   50 continue
1636       end select
1637 !      continue
1638       close (isidep)
1639 !-----------------------------------------------------------------------
1640 ! Calculate the "working" parameters of SC interactions.
1641
1642 !el from module energy - COMMON.INTERACT-------
1643       allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
1644       allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
1645       do i=1,ntyp1
1646         do j=1,ntyp1
1647           aa(i,j)=0.0D0
1648           bb(i,j)=0.0D0
1649           chi(i,j)=0.0D0
1650           sigma(i,j)=0.0D0
1651           r0(i,j)=0.0D0
1652         enddo
1653       enddo
1654 !--------------------------------
1655
1656       do i=2,ntyp
1657         do j=1,i-1
1658           eps(i,j)=eps(j,i)
1659         enddo
1660       enddo
1661       do i=1,ntyp
1662         do j=i,ntyp
1663           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1664           sigma(j,i)=sigma(i,j)
1665           rs0(i,j)=dwa16*sigma(i,j)
1666           rs0(j,i)=rs0(i,j)
1667         enddo
1668       enddo
1669       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
1670        'Working parameters of the SC interactions:',&
1671        '     a    ','     b    ','   augm   ','  sigma ','   r0   ',&
1672        '  chi1   ','   chi2   '
1673       do i=1,ntyp
1674         do j=i,ntyp
1675           epsij=eps(i,j)
1676           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1677             rrij=sigma(i,j)
1678           else
1679             rrij=rr0(i)+rr0(j)
1680           endif
1681           r0(i,j)=rrij
1682           r0(j,i)=rrij
1683           rrij=rrij**expon
1684           epsij=eps(i,j)
1685           sigeps=dsign(1.0D0,epsij)
1686           epsij=dabs(epsij)
1687           aa(i,j)=epsij*rrij*rrij
1688           bb(i,j)=-sigeps*epsij*rrij
1689           aa(j,i)=aa(i,j)
1690           bb(j,i)=bb(i,j)
1691           if (ipot.gt.2) then
1692             sigt1sq=sigma0(i)**2
1693             sigt2sq=sigma0(j)**2
1694             sigii1=sigii(i)
1695             sigii2=sigii(j)
1696             ratsig1=sigt2sq/sigt1sq
1697             ratsig2=1.0D0/ratsig1
1698             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1699             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1700             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1701           else
1702             rsum_max=sigma(i,j)
1703           endif
1704 !         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1705             sigmaii(i,j)=rsum_max
1706             sigmaii(j,i)=rsum_max
1707 !         else
1708 !           sigmaii(i,j)=r0(i,j)
1709 !           sigmaii(j,i)=r0(i,j)
1710 !         endif
1711 !d        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1712           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1713             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1714             augm(i,j)=epsij*r_augm**(2*expon)
1715 !           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1716             augm(j,i)=augm(i,j)
1717           else
1718             augm(i,j)=0.0D0
1719             augm(j,i)=0.0D0
1720           endif
1721           if (lprint) then
1722             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')  &
1723             restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),&
1724             sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1725           endif
1726         enddo
1727       enddo
1728
1729       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
1730       do i=1,ntyp
1731         do j=1,2
1732           bad(i,j)=0.0D0
1733         enddo
1734       enddo
1735 #ifdef CLUSTER
1736 !
1737 ! Define the SC-p interaction constants
1738 !
1739       do i=1,20
1740         do j=1,2
1741           eps_scp(i,j)=-1.5d0
1742           rscp(i,j)=4.0d0
1743         enddo
1744       enddo
1745 #endif
1746
1747 !elwrite(iout,*) "parmread kontrol before oldscp"
1748 !
1749 ! Define the SC-p interaction constants
1750 !
1751 #ifdef OLDSCP
1752       do i=1,20
1753 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1754 ! helix formation)
1755 !       aad(i,1)=0.3D0*4.0D0**12
1756 ! Following line for constants currently implemented
1757 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
1758         aad(i,1)=1.5D0*4.0D0**12
1759 !       aad(i,1)=0.17D0*5.6D0**12
1760         aad(i,2)=aad(i,1)
1761 ! "Soft" SC-p repulsion
1762         bad(i,1)=0.0D0
1763 ! Following line for constants currently implemented
1764 !       aad(i,1)=0.3D0*4.0D0**6
1765 ! "Hard" SC-p repulsion
1766         bad(i,1)=3.0D0*4.0D0**6
1767 !       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1768         bad(i,2)=bad(i,1)
1769 !       aad(i,1)=0.0D0
1770 !       aad(i,2)=0.0D0
1771 !       bad(i,1)=1228.8D0
1772 !       bad(i,2)=1228.8D0
1773       enddo
1774 #else
1775 !
1776 ! 8/9/01 Read the SC-p interaction constants from file
1777 !
1778       do i=1,ntyp
1779         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1780       enddo
1781       do i=1,ntyp
1782         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1783         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1784         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1785         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1786       enddo
1787
1788       if (lprint) then
1789         write (iout,*) "Parameters of SC-p interactions:"
1790         do i=1,20
1791           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
1792            eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1793         enddo
1794       endif
1795 #endif
1796 !
1797 ! Define the constants of the disulfide bridge
1798 !
1799       ebr=-5.50D0
1800 !
1801 ! Old arbitrary potential - commented out.
1802 !
1803 !      dbr= 4.20D0
1804 !      fbr= 3.30D0
1805 !
1806 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1807 ! energy surface of diethyl disulfide.
1808 ! A. Liwo and U. Kozlowska, 11/24/03
1809 !
1810       D0CM = 3.78d0
1811       AKCM = 15.1d0
1812       AKTH = 11.0d0
1813       AKCT = 12.0d0
1814       V1SS =-1.08d0
1815       V2SS = 7.61d0
1816       V3SS = 13.7d0
1817
1818       if (lprint) then
1819       write (iout,'(/a)') "Disulfide bridge parameters:"
1820       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1821       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1822       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1823       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
1824        ' v3ss:',v3ss
1825       endif
1826       return
1827       end subroutine parmread
1828 #ifndef CLUSTER
1829 !-----------------------------------------------------------------------------
1830 ! mygetenv.F
1831 !-----------------------------------------------------------------------------
1832       subroutine mygetenv(string,var)
1833 !
1834 ! Version 1.0
1835 !
1836 ! This subroutine passes the environmental variables to FORTRAN program.
1837 ! If the flags -DMYGETENV and -DMPI are not for compilation, it calls the
1838 ! standard FORTRAN GETENV subroutine. If both flags are set, the subroutine
1839 ! reads the environmental variables from $HOME/.env
1840 !
1841 ! Usage: As for the standard FORTRAN GETENV subroutine.
1842 !
1843 ! Purpose: some versions/installations of MPI do not transfer the environmental
1844 ! variables to slave processors, if these variables are set in the shell script
1845 ! from which mpirun is called.
1846 !
1847 ! A.Liwo, 7/29/01
1848 !
1849 #ifdef MPI
1850       use MPI_data
1851       include "mpif.h"
1852 #endif
1853 !      implicit none
1854       character*(*) :: string,var
1855 #if defined(MYGETENV) && defined(MPI)
1856 !      include "DIMENSIONS.ZSCOPT"
1857 !      include "mpif.h"
1858 !      include "COMMON.MPI"
1859 !el      character*360 ucase
1860 !el      external ucase
1861       character(len=360) :: string1(360),karta
1862       character(len=240) :: home
1863       integer i,n !,ilen
1864 !el      external ilen
1865       call getenv("HOME",home)
1866       open(99,file=home(:ilen(home))//"/.env",status="OLD",err=112)
1867       do while (.true.)
1868         read (99,end=111,err=111,'(a)') karta
1869         do i=1,80
1870           string1(i)=" "
1871         enddo
1872         call split_string(karta,string1,80,n)
1873         if (ucase(string1(1)(:ilen(string1(1)))).eq."SETENV" .and. &
1874          string1(2)(:ilen(string1(2))).eq.string(:ilen(string)) ) then
1875            var=string1(3)
1876            print *,"Processor",me,": ",var(:ilen(var)),&
1877             " assigned to ",string(:ilen(string))
1878            close(99)
1879            return
1880         endif
1881       enddo
1882  111  print *,"Environment variable ",string(:ilen(string))," not set."
1883       close(99)
1884       return
1885  112  print *,"Error opening environment file!"
1886 #else
1887       call getenv(string,var)
1888 #endif
1889       return
1890       end subroutine mygetenv
1891 !-----------------------------------------------------------------------------
1892 ! readrtns.F
1893 !-----------------------------------------------------------------------------
1894       subroutine read_general_data(*)
1895
1896       use control_data, only:indpdb,symetr
1897       use energy_data, only:distchainmax
1898 !      implicit none
1899 !      include "DIMENSIONS"
1900 !      include "DIMENSIONS.ZSCOPT"
1901 !      include "DIMENSIONS.FREE"
1902 !      include "COMMON.TORSION"
1903 !      include "COMMON.INTERACT"
1904 !      include "COMMON.IOUNITS"
1905 !      include "COMMON.TIME1"
1906 !      include "COMMON.PROT"
1907 !      include "COMMON.PROTFILES"
1908 !      include "COMMON.CHAIN"
1909 !      include "COMMON.NAMES"
1910 !      include "COMMON.FFIELD"
1911 !      include "COMMON.ENEPS"
1912 !      include "COMMON.WEIGHTS"
1913 !      include "COMMON.FREE"
1914 !      include "COMMON.CONTROL"
1915 !      include "COMMON.ENERGIES"
1916       character(len=800) :: controlcard
1917       integer :: i,j,k,ii,n_ene_found
1918       integer :: ind,itype1,itype2,itypf,itypsc,itypp
1919 !el      integer ilen
1920 !el      external ilen
1921 !el      character*16 ucase
1922       character(len=16) :: key
1923 !el      external ucase
1924       call card_concat(controlcard,.true.)
1925       call readi(controlcard,"N_ENE",n_eneW,max_eneW)
1926       if (n_eneW.gt.max_eneW) then
1927         write (iout,*) "Error: parameter out of range: N_ENE",n_eneW,&
1928           max_eneW
1929         return 1
1930       endif
1931       call readi(controlcard,"NPARMSET",nparmset,1)
1932 !elwrite(iout,*)"in read_gen data"
1933       separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
1934       call readi(controlcard,"IPARMPRINT",iparmprint,1)
1935       write (iout,*) "PARMPRINT",iparmprint
1936       if (nparmset.gt.max_parm) then
1937         write (iout,*) "Error: parameter out of range: NPARMSET",&
1938           nparmset, Max_Parm
1939         return 1
1940       endif
1941 !elwrite(iout,*)"in read_gen data"
1942       call readi(controlcard,"MAXIT",maxit,5000)
1943       call reada(controlcard,"FIMIN",fimin,1.0d-3)
1944       call readi(controlcard,"ENSEMBLES",ensembles,0)
1945       hamil_rep=index(controlcard,"HAMIL_REP").gt.0
1946       write (iout,*) "Number of energy parameter sets",nparmset
1947       allocate(isampl(nparmset))
1948       call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
1949       write (iout,*) "MaxSlice",MaxSlice
1950       call readi(controlcard,"NSLICE",nslice,1)
1951 !elwrite(iout,*)"in read_gen data"
1952       call flush(iout)
1953       if (nslice.gt.MaxSlice) then
1954         write (iout,*) "Error: parameter out of range: NSLICE",nslice,&
1955           MaxSlice
1956         return 1
1957       endif
1958       write (iout,*) "Frequency of storing conformations",&
1959        (isampl(i),i=1,nparmset)
1960       write (iout,*) "Maxit",maxit," Fimin",fimin
1961       call readi(controlcard,"NQ",nQ,1)
1962       if (nQ.gt.MaxQ) then
1963         write (iout,*) "Error: parameter out of range: NQ",nq,&
1964           maxq
1965         return 1
1966       endif
1967       indpdb=0
1968       if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
1969       call reada(controlcard,"DELTA",delta,1.0d-2)
1970       call readi(controlcard,"EINICHECK",einicheck,2)
1971       call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
1972       call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
1973       call readi(controlcard,"RESCALE",rescale_modeW,1)
1974       check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
1975       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
1976       call readi(controlcard,'SYM',symetr,1)
1977       write (iout,*) "DISTCHAINMAX",distchainmax
1978       write (iout,*) "delta",delta
1979       write (iout,*) "einicheck",einicheck
1980       write (iout,*) "rescale_mode",rescale_modeW
1981       call flush(iout)
1982       bxfile=index(controlcard,"BXFILE").gt.0
1983       cxfile=index(controlcard,"CXFILE").gt.0
1984       if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile) &
1985        bxfile=.true.
1986       histfile=index(controlcard,"HISTFILE").gt.0
1987       histout=index(controlcard,"HISTOUT").gt.0
1988       entfile=index(controlcard,"ENTFILE").gt.0
1989       zscfile=index(controlcard,"ZSCFILE").gt.0
1990       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1991       write (iout,*) "with_dihed_constr ",with_dihed_constr
1992       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
1993       return
1994       end subroutine read_general_data
1995 !------------------------------------------------------------------------------
1996       subroutine read_efree(*)
1997 !
1998 ! Read molecular data
1999 !
2000 !      implicit none
2001 !      include 'DIMENSIONS'
2002 !      include 'DIMENSIONS.ZSCOPT'
2003 !      include 'DIMENSIONS.COMPAR'
2004 !      include 'DIMENSIONS.FREE'
2005 !      include 'COMMON.IOUNITS'
2006 !      include 'COMMON.TIME1'
2007 !      include 'COMMON.SBRIDGE'
2008 !      include 'COMMON.CONTROL'
2009 !      include 'COMMON.CHAIN'
2010 !      include 'COMMON.HEADER'
2011 !      include 'COMMON.GEO'
2012 !      include 'COMMON.FREE'
2013       character(len=320) :: controlcard !,ucase
2014       integer :: iparm,ib,i,j,npars
2015 !el      integer ilen
2016 !el      external ilen
2017
2018       if (hamil_rep) then
2019         npars=1
2020       else
2021         npars=nParmSet
2022       endif
2023
2024 !      call alloc_wham_arrays
2025 !      allocate(nT_h(nParmSet))
2026 !      allocate(replica(nParmSet))
2027 !      allocate(umbrella(nParmSet))
2028 !      allocate(read_iset(nParmSet))
2029 !      allocate(nT_h(nParmSet))
2030
2031       do iparm=1,npars
2032
2033       call card_concat(controlcard,.true.)
2034       call readi(controlcard,'NT',nT_h(iparm),1)
2035       write (iout,*) "iparm",iparm," nt",nT_h(iparm)
2036       call flush(iout)
2037       if (nT_h(iparm).gt.MaxT_h) then
2038         write (iout,*)  "Error: parameter out of range: NT",nT_h(iparm),&
2039           MaxT_h
2040         return 1
2041       endif
2042       replica(iparm)=index(controlcard,"REPLICA").gt.0
2043       umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
2044       read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
2045       write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",&
2046         replica(iparm)," umbrella ",umbrella(iparm),&
2047         " read_iset",read_iset(iparm)
2048       call flush(iout)
2049       do ib=1,nT_h(iparm)
2050         call card_concat(controlcard,.true.)
2051         call readi(controlcard,'NR',nR(ib,iparm),1)
2052         if (umbrella(iparm)) then
2053           nRR(ib,iparm)=1
2054         else
2055           nRR(ib,iparm)=nR(ib,iparm)
2056         endif
2057         if (nR(ib,iparm).gt.MaxR) then
2058           write (iout,*)  "Error: parameter out of range: NR",&
2059             nR(ib,iparm),MaxR
2060         return 1
2061         endif
2062         call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
2063         beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
2064         call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),&
2065           0.0d0)
2066         do i=1,nR(ib,iparm)
2067           call card_concat(controlcard,.true.)
2068           call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,&
2069             100.0d0)
2070           call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,&
2071             0.0d0)
2072         enddo
2073       enddo
2074       do ib=1,nT_h(iparm)
2075         write (iout,*) "ib",ib," beta_h",&
2076           1.0d0/(0.001987*beta_h(ib,iparm))
2077         write (iout,*) "nR",nR(ib,iparm)
2078         write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
2079         do i=1,nR(ib,iparm)
2080           write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),&
2081             "q0",(q0(j,i,ib,iparm),j=1,nQ)
2082         enddo
2083         call flush(iout)
2084       enddo
2085
2086       enddo
2087
2088       if (hamil_rep) then
2089
2090        do iparm=2,nParmSet
2091           nT_h(iparm)=nT_h(1)
2092          do ib=1,nT_h(iparm)
2093            nR(ib,iparm)=nR(ib,1)
2094            if (umbrella(iparm)) then
2095              nRR(ib,iparm)=1
2096            else
2097              nRR(ib,iparm)=nR(ib,1)
2098            endif
2099            beta_h(ib,iparm)=beta_h(ib,1)
2100            do i=1,nR(ib,iparm)
2101              f(i,ib,iparm)=f(i,ib,1)
2102              do j=1,nQ
2103                KH(j,i,ib,iparm)=KH(j,i,ib,1)
2104                Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
2105              enddo
2106            enddo
2107            replica(iparm)=replica(1)
2108            umbrella(iparm)=umbrella(1)
2109            read_iset(iparm)=read_iset(1)
2110          enddo
2111        enddo
2112
2113       endif
2114
2115       return
2116       end subroutine read_efree
2117 !-----------------------------------------------------------------------------
2118       subroutine read_protein_data(*)
2119 !      implicit none
2120 !      include "DIMENSIONS"
2121 !      include "DIMENSIONS.ZSCOPT"
2122 !      include "DIMENSIONS.FREE"
2123 #ifdef MPI
2124       use MPI_data
2125       include "mpif.h"
2126       integer :: IERROR,ERRCODE!,STATUS(MPI_STATUS_SIZE)
2127 !      include "COMMON.MPI"
2128 #endif
2129 !      include "COMMON.CHAIN"
2130 !      include "COMMON.IOUNITS"
2131 !      include "COMMON.PROT"
2132 !      include "COMMON.PROTFILES"
2133 !      include "COMMON.NAMES"
2134 !      include "COMMON.FREE"
2135 !      include "COMMON.OBCINKA"
2136       character(len=64) :: nazwa
2137       character(len=16000) :: controlcard
2138       integer :: i,ii,ib,iR,iparm,nthr,npars !,ilen,iroof
2139 !el      external ilen,iroof
2140       if (hamil_rep) then
2141         npars=1
2142       else
2143         npars=nparmset
2144       endif
2145
2146       do iparm=1,npars
2147
2148 ! Read names of files with conformation data.
2149       if (replica(iparm)) then
2150         nthr = 1
2151       else
2152         nthr = nT_h(iparm)
2153       endif
2154       do ib=1,nthr
2155       do ii=1,nRR(ib,iparm)
2156       write (iout,*) "Parameter set",iparm," temperature",ib,&
2157        " window",ii
2158       call flush(iout)
2159       call card_concat(controlcard,.true.)
2160       write (iout,*) controlcard(:ilen(controlcard))
2161       call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
2162       call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
2163       call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
2164       call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
2165       call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),&
2166        maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
2167       call reada(controlcard,"TIME_START",&
2168         time_start_collect(ii,ib,iparm),0.0d0)
2169       call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),&
2170         1.0d10)
2171       write (iout,*) "rec_start",rec_start(ii,ib,iparm),&
2172        " rec_end",rec_end(ii,ib,iparm)
2173       write (iout,*) "time_start",time_start_collect(ii,ib,iparm),&
2174        " time_end",time_end_collect(ii,ib,iparm)
2175       call flush(iout)
2176       if (replica(iparm)) then
2177         call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
2178         write (iout,*) "Number of trajectories",totraj(ii,iparm)
2179         call flush(iout)
2180       endif
2181       if (nfile_bin(ii,ib,iparm).lt.2 &
2182           .and. nfile_asc(ii,ib,iparm).eq.0 &
2183           .and. nfile_cx(ii,ib,iparm).eq.0) then
2184         write (iout,*) "Error - no action specified!"
2185         return 1
2186       endif
2187       if (nfile_bin(ii,ib,iparm).gt.0) then
2188         call card_concat(controlcard,.false.)
2189         call split_string(controlcard,protfiles(1,1,ii,ib,iparm),&
2190          maxfile_prot,nfile_bin(ii,ib,iparm))
2191 #ifdef DEBUG
2192         write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
2193         write(iout,*) (protfiles(i,1,ii,ib,iparm),&
2194           i=1,nfile_bin(ii,ib,iparm))
2195 #endif
2196       endif
2197       if (nfile_asc(ii,ib,iparm).gt.0) then
2198         call card_concat(controlcard,.false.)
2199         call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
2200          maxfile_prot,nfile_asc(ii,ib,iparm))
2201 #ifdef DEBUG
2202         write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
2203         write(iout,*) (protfiles(i,2,ii,ib,iparm),&
2204           i=1,nfile_asc(ii,ib,iparm))
2205 #endif
2206       else if (nfile_cx(ii,ib,iparm).gt.0) then
2207         call card_concat(controlcard,.false.)
2208         call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
2209          maxfile_prot,nfile_cx(ii,ib,iparm))
2210 #ifdef DEBUG
2211         write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
2212         write(iout,*) (protfiles(i,2,ii,ib,iparm),&
2213          i=1,nfile_cx(ii,ib,iparm))
2214 #endif
2215       endif
2216       call flush(iout)
2217       enddo
2218       enddo
2219
2220       enddo
2221       return
2222       end subroutine read_protein_data
2223 !-------------------------------------------------------------------------------
2224       subroutine readsss(rekord,lancuch,wartosc,default)
2225 !      implicit none
2226       character*(*) :: rekord,lancuch,wartosc,default
2227       character(len=80) :: aux
2228       integer :: lenlan,lenrec,iread,ireade
2229 !el      external ilen
2230 !el      logical iblnk
2231 !el      external iblnk
2232       lenlan=ilen(lancuch)
2233       lenrec=ilen(rekord)
2234       iread=index(rekord,lancuch(:lenlan)//"=")
2235 !      print *,"rekord",rekord," lancuch",lancuch
2236 !      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
2237       if (iread.eq.0) then
2238         wartosc=default
2239         return
2240       endif
2241       iread=iread+lenlan+1
2242 !      print *,"iread",iread
2243 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2244       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
2245         iread=iread+1
2246 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2247       enddo
2248 !      print *,"iread",iread
2249       if (iread.gt.lenrec) then
2250          wartosc=default
2251         return
2252       endif
2253       ireade=iread+1
2254 !      print *,"ireade",ireade
2255       do while (ireade.lt.lenrec .and. &
2256          .not.iblnk(rekord(ireade:ireade)))
2257         ireade=ireade+1
2258       enddo
2259       wartosc=rekord(iread:ireade)
2260       return
2261       end subroutine readsss
2262 !----------------------------------------------------------------------------
2263       subroutine multreads(rekord,lancuch,tablica,dim,default)
2264 !      implicit none
2265       integer :: dim,i
2266       character*(*) rekord,lancuch,tablica(dim),default
2267       character(len=80) :: aux
2268       integer :: lenlan,lenrec,iread,ireade
2269 !el      external ilen
2270 !el      logical iblnk
2271 !el      external iblnk
2272       do i=1,dim
2273         tablica(i)=default
2274       enddo
2275       lenlan=ilen(lancuch)
2276       lenrec=ilen(rekord)
2277       iread=index(rekord,lancuch(:lenlan)//"=")
2278 !      print *,"rekord",rekord," lancuch",lancuch
2279 !      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
2280       if (iread.eq.0) return
2281       iread=iread+lenlan+1
2282       do i=1,dim
2283 !      print *,"iread",iread
2284 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2285       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
2286         iread=iread+1
2287 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2288       enddo
2289 !      print *,"iread",iread
2290       if (iread.gt.lenrec) return
2291       ireade=iread+1
2292 !      print *,"ireade",ireade
2293       do while (ireade.lt.lenrec .and. &
2294          .not.iblnk(rekord(ireade:ireade)))
2295         ireade=ireade+1
2296       enddo
2297       tablica(i)=rekord(iread:ireade)
2298       iread=ireade+1
2299       enddo
2300       end subroutine multreads
2301 !----------------------------------------------------------------------------
2302       subroutine split_string(rekord,tablica,dim,nsub)
2303 !      implicit none
2304       integer :: dim,nsub,i,ii,ll,kk
2305       character*(*) tablica(dim)
2306       character*(*) rekord
2307 !el      integer ilen
2308 !el      external ilen
2309       do i=1,dim
2310         tablica(i)=" "
2311       enddo
2312       ii=1
2313       ll = ilen(rekord)
2314       nsub=0
2315       do i=1,dim
2316 ! Find the start of term name
2317         kk = 0
2318         do while (ii.le.ll .and. rekord(ii:ii).eq." ")
2319           ii = ii+1
2320         enddo
2321 ! Parse the name into TABLICA(i) until blank found
2322         do while (ii.le.ll .and. rekord(ii:ii).ne." ")
2323           kk = kk+1
2324           tablica(i)(kk:kk)=rekord(ii:ii)
2325           ii = ii+1
2326         enddo
2327         if (kk.gt.0) nsub=nsub+1
2328         if (ii.gt.ll) return
2329       enddo
2330       return
2331       end subroutine split_string
2332 !--------------------------------------------------------------------------------
2333 ! readrtns_compar.F
2334 !--------------------------------------------------------------------------------
2335       subroutine read_compar
2336 !
2337 ! Read molecular data
2338 !
2339       use conform_compar, only:alloc_compar_arrays
2340       use control_data, only:pdbref
2341       use geometry_data, only:deg2rad,rad2deg
2342 !      implicit none
2343 !      include 'DIMENSIONS'
2344 !      include 'DIMENSIONS.ZSCOPT'
2345 !      include 'DIMENSIONS.COMPAR'
2346 !      include 'DIMENSIONS.FREE'
2347 !      include 'COMMON.IOUNITS'
2348 !      include 'COMMON.TIME1'
2349 !      include 'COMMON.SBRIDGE'
2350 !      include 'COMMON.CONTROL'
2351 !      include 'COMMON.COMPAR'
2352 !      include 'COMMON.CHAIN'
2353 !      include 'COMMON.HEADER'
2354 !      include 'COMMON.GEO'
2355 !      include 'COMMON.FREE'
2356       character(len=320) :: controlcard !,ucase
2357       character(len=64) :: wfile
2358 !el      integer ilen
2359 !el      external ilen
2360       integer :: i,j,k
2361 !elwrite(iout,*)"jestesmy w read_compar"
2362       call card_concat(controlcard,.true.)
2363       pdbref=(index(controlcard,'PDBREF').gt.0)
2364       call reada(controlcard,'CUTOFF_UP',rmscut_base_up,4.0d0)
2365       call reada(controlcard,'CUTOFF_LOW',rmscut_base_low,3.0d0)
2366       call reada(controlcard,'RMSUP_LIM',rmsup_lim,4.0d0)
2367       call reada(controlcard,'RMSUPUP_LIM',rmsupup_lim,7.5d0)
2368       verbose = index(controlcard,"VERBOSE").gt.0
2369       lgrp=index(controlcard,"STATIN").gt.0
2370       lgrp_out=index(controlcard,"STATOUT").gt.0
2371       merge_helices=index(controlcard,"DONT_MERGE_HELICES").eq.0
2372       binary = index(controlcard,"BINARY").gt.0
2373       rmscut_base_up=rmscut_base_up/50
2374       rmscut_base_low=rmscut_base_low/50
2375       call reada(controlcard,"FRAC_SEC",frac_sec,0.66666666d0)
2376       call readi(controlcard,'NLEVEL',nlevel,1)
2377       if (nlevel.lt.0) then
2378         allocate(nfrag(2))
2379         call alloc_compar_arrays(maxfrag,1)
2380         goto 121
2381       else
2382         allocate(nfrag(nlevel))
2383       endif
2384 ! Read the data pertaining to elementary fragments (level 1)
2385       call readi(controlcard,'NFRAG',nfrag(1),0)
2386       write(iout,*)"nfrag(1)",nfrag(1)
2387       call alloc_compar_arrays(nfrag(1),nlevel)
2388       do j=1,nfrag(1)
2389         call card_concat(controlcard,.true.)
2390         write (iout,*) controlcard(:ilen(controlcard))
2391         call readi(controlcard,'NPIECE',npiece(j,1),0)
2392         call readi(controlcard,'N_SHIFT1',n_shift(1,j,1),0)
2393         call readi(controlcard,'N_SHIFT2',n_shift(2,j,1),0)
2394         call reada(controlcard,'ANGCUT',ang_cut(j),50.0d0)
2395         call reada(controlcard,'MAXANG',ang_cut1(j),360.0d0)
2396         call reada(controlcard,'FRAC_MIN',frac_min(j),0.666666d0)
2397         call reada(controlcard,'NC_FRAC',nc_fragm(j,1),0.5d0)
2398         call readi(controlcard,'NC_REQ',nc_req_setf(j,1),0)
2399         call readi(controlcard,'RMS',irms(j,1),0)
2400         call readi(controlcard,'LOCAL',iloc(j),1)
2401         call readi(controlcard,'ELCONT',ielecont(j,1),1)
2402         if (ielecont(j,1).eq.0) then
2403           call readi(controlcard,'SCCONT',isccont(j,1),1)
2404         endif
2405         ang_cut(j)=ang_cut(j)*deg2rad
2406         ang_cut1(j)=ang_cut1(j)*deg2rad
2407         do k=1,npiece(j,1)
2408           call card_concat(controlcard,.true.)
2409           call readi(controlcard,'IFRAG1',ifrag(1,k,j),0)
2410           call readi(controlcard,'IFRAG2',ifrag(2,k,j),0)
2411         enddo
2412         write(iout,*)"j",j," npiece",npiece(j,1)," ifrag",&
2413           (ifrag(1,k,j),ifrag(2,k,j),&
2414          k=1,npiece(j,1))," ang_cut",ang_cut(j)*rad2deg,&
2415           " ang_cut1",ang_cut1(j)*rad2deg
2416         write(iout,*)"n_shift",n_shift(1,j,1),n_shift(2,j,1)
2417         write(iout,*)"nc_frac",nc_fragm(j,1)," nc_req",nc_req_setf(j,1)
2418         write(iout,*)"irms",irms(j,1)," ielecont",ielecont(j,1),&
2419           " ilocal",iloc(j)," isccont",isccont(j,1)
2420       enddo
2421 ! Read data pertaning to higher levels
2422       do i=2,nlevel
2423         call card_concat(controlcard,.true.)
2424         call readi(controlcard,'NFRAG',NFRAG(i),0)
2425         write (iout,*) "i",i," nfrag",nfrag(i)
2426         do j=1,nfrag(i)
2427           call card_concat(controlcard,.true.)
2428           if (i.eq.2) then
2429             call readi(controlcard,'ELCONT',ielecont(j,i),0)
2430             if (ielecont(j,i).eq.0) then
2431               call readi(controlcard,'SCCONT',isccont(j,i),1)
2432             endif
2433             call readi(controlcard,'RMS',irms(j,i),0)
2434           else
2435             ielecont(j,i)=0
2436             isccont(j,i)=0
2437             irms(j,i)=1
2438           endif
2439           call readi(controlcard,'NPIECE',npiece(j,i),0)
2440           call readi(controlcard,'N_SHIFT1',n_shift(1,j,i),0)
2441           call readi(controlcard,'N_SHIFT2',n_shift(2,j,i),0)
2442           call multreadi(controlcard,'IPIECE',ipiece(1,j,i),&
2443             npiece(j,i),0)
2444           call reada(controlcard,'NC_FRAC',nc_fragm(j,i),0.5d0)
2445           call readi(controlcard,'NC_REQ',nc_req_setf(j,i),0)
2446           write(iout,*) "j",j," npiece",npiece(j,i)," n_shift",&
2447             n_shift(1,j,i),n_shift(2,j,i)," ielecont",ielecont(j,i),&
2448             " isccont",isccont(j,i)," irms",irms(j,i)
2449           write(iout,*) "ipiece",(ipiece(k,j,i),k=1,npiece(j,i))
2450           write(iout,*)"n_shift",n_shift(1,j,i),n_shift(2,j,i)
2451           write(iout,*)"nc_frac",nc_fragm(j,i),&
2452            " nc_req",nc_req_setf(j,i)
2453         enddo
2454       enddo
2455       if (binary) write (iout,*) "Classes written in binary format."
2456       return
2457   121 continue
2458       call reada(controlcard,'ANGCUT_HEL',angcut_hel,50.0d0)
2459       call reada(controlcard,'MAXANG_HEL',angcut1_hel,60.0d0)
2460       call reada(controlcard,'ANGCUT_BET',angcut_bet,90.0d0)
2461       call reada(controlcard,'MAXANG_BET',angcut1_bet,360.0d0)
2462       call reada(controlcard,'ANGCUT_STRAND',angcut_strand,90.0d0)
2463       call reada(controlcard,'MAXANG_STRAND',angcut1_strand,60.0d0)
2464       call reada(controlcard,'FRAC_MIN',frac_min_set,0.666666d0)
2465       call reada(controlcard,'NC_FRAC_HEL',ncfrac_hel,0.5d0)
2466       call readi(controlcard,'NC_REQ_HEL',ncreq_hel,0)
2467       call reada(controlcard,'NC_FRAC_BET',ncfrac_bet,0.5d0)
2468       call reada(controlcard,'NC_FRAC_PAIR',ncfrac_pair,0.3d0)
2469       call readi(controlcard,'NC_REQ_BET',ncreq_bet,0)
2470       call readi(controlcard,'NC_REQ_PAIR',ncreq_pair,0)
2471       call readi(controlcard,'NSHIFT_HEL',nshift_hel,3)
2472       call readi(controlcard,'NSHIFT_BET',nshift_bet,3)
2473       call readi(controlcard,'NSHIFT_STRAND',nshift_strand,3)
2474       call readi(controlcard,'NSHIFT_PAIR',nshift_pair,3)
2475       call readi(controlcard,'RMS_SINGLE',irms_single,0)
2476       call readi(controlcard,'CONT_SINGLE',icont_single,1)
2477       call readi(controlcard,'LOCAL_SINGLE',iloc_single,1)
2478       call readi(controlcard,'RMS_PAIR',irms_pair,0)
2479       call readi(controlcard,'CONT_PAIR',icont_pair,1)
2480       call readi(controlcard,'SPLIT_BET',isplit_bet,0)
2481       angcut_hel=angcut_hel*deg2rad
2482       angcut1_hel=angcut1_hel*deg2rad
2483       angcut_bet=angcut_bet*deg2rad
2484       angcut1_bet=angcut1_bet*deg2rad
2485       angcut_strand=angcut_strand*deg2rad
2486       angcut1_strand=angcut1_strand*deg2rad
2487       write (iout,*) "Automatic detection of structural elements"
2488       write (iout,*) 'NC_FRAC_HEL',ncfrac_hel,' NC_REQ_HEL',ncreq_hel,&
2489                      ' NC_FRAC_BET',ncfrac_bet,' NC_REQ_BET',ncreq_bet,&
2490                  ' RMS_SINGLE',irms_single,' CONT_SINGLE',icont_single,&
2491                  ' NC_FRAC_PAIR',ncfrac_pair,' NC_REQ_PAIR',ncreq_pair,&
2492         ' RMS_PAIR',irms_pair,' CONT_PAIR',icont_pair,&
2493         ' SPLIT_BET',isplit_bet
2494       write (iout,*) 'NSHIFT_HEL',nshift_hel,' NSHIFT_BET',nshift_bet,&
2495         ' NSHIFT_STRAND',nshift_strand,' NSHIFT_PAIR',nshift_pair
2496       write (iout,*) 'ANGCUT_HEL',angcut_hel*rad2deg,&
2497         ' MAXANG_HEL',angcut1_hel*rad2deg
2498       write (iout,*) 'ANGCUT_BET',angcut_bet*rad2deg,&
2499                      ' MAXANG_BET',angcut1_bet*rad2deg
2500       write (iout,*) 'ANGCUT_STRAND',angcut_strand*rad2deg,&
2501                      ' MAXANG_STRAND',angcut1_strand*rad2deg
2502       write (iout,*) 'FRAC_MIN',frac_min_set
2503       return
2504       end subroutine read_compar
2505 !--------------------------------------------------------------------------------
2506 ! read_ref_str.F
2507 !--------------------------------------------------------------------------------
2508       subroutine read_ref_structure(*)
2509 !
2510 ! Read the reference structure from the PDB file or from a PDB file or in the form of the dihedral
2511 ! angles.
2512 !
2513       use control_data, only:pdbref
2514       use geometry_data, only:nres,cref,c,dc,nsup,dc_norm,nend_sup,&
2515                               nstart_sup,nstart_seq,nperm,nres0
2516       use energy_data, only:nct,nnt,icont_ref,ncont_ref,itype
2517       use compare, only:seq_comp !,contact,elecont
2518       use geometry, only:chainbuild,dist
2519       use io_config, only:readpdb
2520 !
2521       use conform_compar, only:contact,elecont
2522 !      implicit none
2523 !      include 'DIMENSIONS'
2524 !      include 'DIMENSIONS.ZSCOPT'
2525 !      include 'DIMENSIONS.COMPAR'
2526 !      include 'COMMON.IOUNITS'
2527 !      include 'COMMON.GEO'
2528 !      include 'COMMON.VAR'
2529 !      include 'COMMON.INTERACT'
2530 !      include 'COMMON.LOCAL'
2531 !      include 'COMMON.NAMES'
2532 !      include 'COMMON.CHAIN'
2533 !      include 'COMMON.FFIELD'
2534 !      include 'COMMON.SBRIDGE'
2535 !      include 'COMMON.HEADER'
2536 !      include 'COMMON.CONTROL'
2537 !      include 'COMMON.CONTACTS1'
2538 !      include 'COMMON.PEPTCONT'
2539 !      include 'COMMON.TIME1'
2540 !      include 'COMMON.COMPAR'
2541       character(len=4) :: sequence(nres)
2542 !el      integer rescode
2543 !el      real(kind=8) :: x(maxvar)
2544       integer :: itype_pdb(nres)
2545 !el      logical seq_comp
2546       integer :: i,j,k,nres_pdb,iaux
2547       real(kind=8) :: ddsc !el,dist
2548       integer :: kkk !,ilen
2549 !el      external ilen
2550 !
2551       nres0=nres
2552       write (iout,*) "pdbref",pdbref
2553       if (pdbref) then
2554         read(inp,'(a)') pdbfile
2555         write (iout,'(2a,1h.)') 'PDB data will be read from file ',&
2556           pdbfile(:ilen(pdbfile))
2557         open(ipdbin,file=pdbfile,status='old',err=33)
2558         goto 34
2559   33    write (iout,'(a)') 'Error opening PDB file.'
2560         return 1
2561   34    continue
2562         do i=1,nres
2563           itype_pdb(i)=itype(i)
2564         enddo
2565
2566         call readpdb
2567
2568         do i=1,nres
2569           iaux=itype_pdb(i)
2570           itype_pdb(i)=itype(i)
2571           itype(i)=iaux
2572         enddo
2573         close (ipdbin)
2574         do kkk=1,nperm
2575         nres_pdb=nres
2576         nres=nres0
2577         nstart_seq=nnt
2578         if (nsup.le.(nct-nnt+1)) then
2579           do i=0,nct-nnt+1-nsup
2580             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),&
2581               nsup)) then
2582               do j=nnt+nsup-1,nnt,-1
2583                 do k=1,3
2584                   cref(k,nres+j+i,kkk)=cref(k,nres_pdb+j,kkk)
2585                 enddo
2586               enddo
2587               do j=nnt+nsup-1,nnt,-1
2588                 do k=1,3
2589                   cref(k,j+i,kkk)=cref(k,j,kkk)
2590                 enddo
2591                 phi_ref(j+i)=phi_ref(j)
2592                 theta_ref(j+i)=theta_ref(j)
2593                 alph_ref(j+i)=alph_ref(j)
2594                 omeg_ref(j+i)=omeg_ref(j)
2595               enddo
2596 #ifdef DEBUG
2597               do j=nnt,nct
2598                 write (iout,'(i5,3f10.5,5x,3f10.5)') &
2599                   j,(cref(k,j,kkk),k=1,3),(cref(k,j+nres,kkk),k=1,3)
2600               enddo
2601 #endif
2602               nstart_seq=nnt+i
2603               nstart_sup=nnt+i
2604               goto 111
2605             endif
2606           enddo
2607           write (iout,'(a)') &
2608                   'Error - sequences to be superposed do not match.'
2609           return 1
2610         else
2611           do i=0,nsup-(nct-nnt+1)
2612             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),&
2613               nct-nnt+1)) &
2614             then
2615               nstart_sup=nstart_sup+i
2616               nsup=nct-nnt+1
2617               goto 111
2618             endif
2619           enddo
2620           write (iout,'(a)') &
2621                   'Error - sequences to be superposed do not match.'
2622         endif
2623         enddo
2624   111   continue
2625         write (iout,'(a,i5)') &
2626          'Experimental structure begins at residue',nstart_seq
2627       else
2628         call read_angles(inp,*38)
2629         goto 39
2630    38   write (iout,'(a)') 'Error reading reference structure.'
2631         return 1
2632    39   call chainbuild
2633         kkk=1
2634         nstart_sup=nnt
2635         nstart_seq=nnt
2636         nsup=nct-nnt+1
2637         do i=1,2*nres
2638           do j=1,3
2639             cref(j,i,kkk)=c(j,i)
2640           enddo
2641         enddo
2642       endif
2643       nend_sup=nstart_sup+nsup-1
2644       do i=1,2*nres
2645         do j=1,3
2646           c(j,i)=cref(j,i,kkk)
2647         enddo
2648       enddo
2649       do i=1,nres
2650         do j=1,3
2651           dc(j,nres+i)=cref(j,nres+i,kkk)-cref(j,i,kkk)
2652         enddo
2653         if (itype(i).ne.10) then
2654           ddsc = dist(i,nres+i)
2655           do j=1,3
2656             dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
2657           enddo
2658         else
2659           do j=1,3
2660             dc_norm(j,nres+i)=0.0d0
2661           enddo
2662         endif
2663 !        write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
2664 !         " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
2665 !         dc_norm(3,nres+i)**2
2666         do j=1,3
2667           dc(j,i)=c(j,i+1)-c(j,i)
2668         enddo
2669         ddsc = dist(i,i+1)
2670         do j=1,3
2671           dc_norm(j,i)=dc(j,i)/ddsc
2672         enddo
2673       enddo
2674 !      print *,"Calling contact"
2675       call contact(.true.,ncont_ref,icont_ref(1,1),&
2676         nstart_sup,nend_sup)
2677 !      print *,"Calling elecont"
2678       call elecont(.true.,ncont_pept_ref,&
2679          icont_pept_ref(1,1),&
2680          nstart_sup,nend_sup)
2681        write (iout,'(a,i3,a,i3,a,i3,a)') &
2682           'Number of residues to be superposed:',nsup,&
2683           ' (from residue',nstart_sup,' to residue',&
2684           nend_sup,').'
2685       return
2686       end subroutine read_ref_structure
2687 !--------------------------------------------------------------------------------
2688 ! geomout.F
2689 !--------------------------------------------------------------------------------
2690       subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev)
2691
2692       use geometry_data, only:nres,c
2693       use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype
2694 !      implicit real*8 (a-h,o-z)
2695 !      include 'DIMENSIONS'
2696 !      include 'DIMENSIONS.ZSCOPT'
2697 !      include 'COMMON.CHAIN'
2698 !      include 'COMMON.INTERACT'
2699 !      include 'COMMON.NAMES'
2700 !      include 'COMMON.IOUNITS'
2701 !      include 'COMMON.HEADER'
2702 !      include 'COMMON.SBRIDGE'
2703       character(len=50) :: tytul
2704       character(len=1),dimension(10) :: chainid=reshape((/'A','B','C',&
2705                       'D','E','F','G','H','I','J'/),shape(chainid))
2706       integer,dimension(nres) :: ica !(maxres)
2707       real(kind=8) :: temp,efree,etot,entropy,rmsdev
2708       integer :: ii,i,j,iti,ires,iatom,ichain
2709       write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')&
2710         ii,temp,rmsdev
2711       write (ipdb,'("REMARK DIMENSIONLESS FREE ENERGY",1pe15.5)') &
2712         efree
2713       write (ipdb,'("REMARK ENERGY",1pe15.5," ENTROPY",1pe15.5)') &
2714         etot,entropy
2715       iatom=0
2716       ichain=1
2717       ires=0
2718       do i=nnt,nct
2719         iti=itype(i)
2720         if (iti.eq.ntyp1) then
2721           ichain=ichain+1
2722           ires=0
2723           write (ipdb,'(a)') 'TER'
2724         else
2725         ires=ires+1
2726         iatom=iatom+1
2727         ica(i)=iatom
2728         write (ipdb,10) iatom,restyp(iti),chainid(ichain),&
2729            ires,(c(j,i),j=1,3)
2730         if (iti.ne.10) then
2731           iatom=iatom+1
2732           write (ipdb,20) iatom,restyp(iti),chainid(ichain),&
2733             ires,(c(j,nres+i),j=1,3)
2734         endif
2735         endif
2736       enddo
2737       write (ipdb,'(a)') 'TER'
2738       do i=nnt,nct-1
2739         if (itype(i).eq.ntyp1) cycle
2740         if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
2741           write (ipdb,30) ica(i),ica(i+1)
2742         else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
2743           write (ipdb,30) ica(i),ica(i+1),ica(i)+1
2744         else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
2745           write (ipdb,30) ica(i),ica(i)+1
2746         endif
2747       enddo
2748       if (itype(nct).ne.10) then
2749         write (ipdb,30) ica(nct),ica(nct)+1
2750       endif
2751       do i=1,nss
2752         write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
2753       enddo
2754       write (ipdb,'(a6)') 'ENDMDL'
2755   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
2756   20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
2757   30  FORMAT ('CONECT',8I5)
2758       return
2759       end subroutine pdboutW
2760 #endif
2761 !------------------------------------------------------------------------------
2762       end module io_wham
2763 !-----------------------------------------------------------------------------
2764 !-----------------------------------------------------------------------------
2765