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: maxterm,maxlor,maxterm_sccor,& !maxtor
 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(-nthetyp-1:nthetyp+1,&
 687         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 688 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 689       allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
 690         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 691 !(maxtheterm,-maxthetyp1:maxthetyp1,&
 692 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 693       allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 694         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 695       allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 696         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 697       allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 698         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 699       allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 700         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 701 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
 702 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 703       allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
 704         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 705       allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
 706         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,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       aa0thet(:,:,:,:)=0.0d0
 717       aathet(:,:,:,:,:)=0.0d0
 718       bbthet(:,:,:,:,:,:)=0.0d0
 719       ccthet(:,:,:,:,:,:)=0.0d0
 720       ddthet(:,:,:,:,:,:)=0.0d0
 721       eethet(:,:,:,:,:,:)=0.0d0
 722       ffthet(:,:,:,:,:,:,:)=0.0d0
 723       ggthet(:,:,:,:,:,:,:)=0.0d0
 724
 725       do iblock=1,2
 726       do i=0,nthetyp
 727         do j=-nthetyp,nthetyp
 728           do k=-nthetyp,nthetyp
 729             read (ithep,'(6a)') res1
 730             read (ithep,*) aa0thet(i,j,k,iblock)
 731             read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
 732             read (ithep,*) &
 733              ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
 734               (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
 735               (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
 736               (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
 737               ll=1,ntheterm2)
 738             read (ithep,*) &
 739             (((ffthet(llll,lll,ll,i,j,k,iblock),&
 740                ffthet(lll,llll,ll,i,j,k,iblock),&
 741                ggthet(llll,lll,ll,i,j,k,iblock),&
 742                ggthet(lll,llll,ll,i,j,k,iblock),&
 743                llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 744           enddo
 745         enddo
 746       enddo
 747 !
 748 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
 749 ! coefficients of theta-and-gamma-dependent terms are zero.
 750 !
 751       do i=1,nthetyp
 752         do j=1,nthetyp
 753           do l=1,ntheterm
 754             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
 755             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
 756           enddo
 757           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
 758           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
 759         enddo
 760         do l=1,ntheterm
 761           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 762         enddo
 763         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 764       enddo
 765       enddo
 766 ! Substitution for D aminoacids from symmetry.
 767       do iblock=1,2
 768       do i=-nthetyp,0
 769         do j=-nthetyp,nthetyp
 770           do k=-nthetyp,nthetyp
 771            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
 772            do l=1,ntheterm
 773            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
 774            enddo
 775            do ll=1,ntheterm2
 776             do lll=1,nsingle
 777             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
 778             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
 779             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
 780             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
 781             enddo
 782           enddo
 783           do ll=1,ntheterm3
 784            do lll=2,ndouble
 785             do llll=1,lll-1
 786             ffthet(llll,lll,ll,i,j,k,iblock)= &
 787             ffthet(llll,lll,ll,-i,-j,-k,iblock)
 788             ffthet(lll,llll,ll,i,j,k,iblock)= &
 789             ffthet(lll,llll,ll,-i,-j,-k,iblock)
 790             ggthet(llll,lll,ll,i,j,k,iblock)= &
 791             -ggthet(llll,lll,ll,-i,-j,-k,iblock)
 792             ggthet(lll,llll,ll,i,j,k,iblock)= &
 793             -ggthet(lll,llll,ll,-i,-j,-k,iblock)
 794             enddo !ll
 795            enddo  !lll
 796           enddo   !llll
 797          enddo    !k
 798         enddo     !j
 799        enddo      !i
 800       enddo       !iblock
 801
 802 !
 803 ! Control printout of the coefficients of virtual-bond-angle potentials
 804 !
 805 do iblock=1,2
 806       if (lprint) then
 807         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 808         do i=1,nthetyp+1
 809           do j=1,nthetyp+1
 810             do k=1,nthetyp+1
 811               write (iout,'(//4a)') &
 812                'Type ',onelett(i),onelett(j),onelett(k)
 813               write (iout,'(//a,10x,a)') " l","a[l]"
 814               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
 815               write (iout,'(i2,1pe15.5)') &
 816                  (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
 817             do l=1,ntheterm2
 818               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))') &
 819                 "b",l,"c",l,"d",l,"e",l
 820               do m=1,nsingle
 821                 write (iout,'(i2,4(1pe15.5))') m,&
 822                 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
 823                 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
 824               enddo
 825             enddo
 826             do l=1,ntheterm3
 827               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))') &
 828                 "f+",l,"f-",l,"g+",l,"g-",l
 829               do m=2,ndouble
 830                 do n=1,m-1
 831                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
 832                     ffthet(n,m,l,i,j,k,iblock),&
 833                     ffthet(m,n,l,i,j,k,iblock),&
 834                     ggthet(n,m,l,i,j,k,iblock),&
 835                     ggthet(m,n,l,i,j,k,iblock)
 836                 enddo
 837               enddo
 838             enddo
 839           enddo
 840         enddo
 841       enddo
 842       call flush(iout)
 843       endif
 844 enddo
 845 #endif
 846 !-------------------------------------------
 847       allocate(nlob(ntyp1)) !(ntyp1)
 848       allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
 849       allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
 850       allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
 851
 852       do i=1,ntyp
 853         do j=1,maxlob
 854           bsc(j,i)=0.0D0
 855           nlob(i)=0
 856         enddo
 857       enddo
 858       nlob(ntyp1)=0
 859       dsc(ntyp1)=0.0D0
 860
 861       do i=-ntyp,ntyp
 862         do j=1,maxlob
 863           do k=1,3
 864             censc(k,j,i)=0.0D0
 865           enddo
 866           do k=1,3
 867             do l=1,3
 868               gaussc(l,k,j,i)=0.0D0
 869             enddo
 870           enddo
 871         enddo
 872       enddo
 873
 874 #ifdef CRYST_SC
 875 !
 876 ! Read the parameters of the probability distribution/energy expression
 877 ! of the side chains.
 878 !
 879       do i=1,ntyp
 880 !c      write (iout,*) "tu dochodze",i
 881         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
 882         if (i.eq.10) then
 883           dsc_inv(i)=0.0D0
 884         else
 885           dsc_inv(i)=1.0D0/dsc(i)
 886         endif
 887         if (i.ne.10) then
 888         do j=1,nlob(i)
 889           do k=1,3
 890             do l=1,3
 891               blower(l,k,j)=0.0D0
 892             enddo
 893           enddo
 894         enddo
 895         bsc(1,i)=0.0D0
 896         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
 897         censc(1,1,-i)=censc(1,1,i)
 898         censc(2,1,-i)=censc(2,1,i)
 899         censc(3,1,-i)=-censc(3,1,i)
 900         do j=2,nlob(i)
 901           read (irotam,*) bsc(j,i)
 902           read (irotam,*) (censc(k,j,i),k=1,3),&
 903                                        ((blower(k,l,j),l=1,k),k=1,3)
 904         censc(1,j,-i)=censc(1,j,i)
 905         censc(2,j,-i)=censc(2,j,i)
 906         censc(3,j,-i)=-censc(3,j,i)
 907 ! BSC is amplitude of Gaussian
 908         enddo
 909         do j=1,nlob(i)
 910           do k=1,3
 911             do l=1,k
 912               akl=0.0D0
 913               do m=1,3
 914                 akl=akl+blower(k,m,j)*blower(l,m,j)
 915               enddo
 916               gaussc(k,l,j,i)=akl
 917               gaussc(l,k,j,i)=akl
 918              if (((k.eq.3).and.(l.ne.3)) &
 919               .or.((l.eq.3).and.(k.ne.3))) then
 920                 gaussc(k,l,j,-i)=-akl
 921                 gaussc(l,k,j,-i)=-akl
 922               else
 923                 gaussc(k,l,j,-i)=akl
 924                 gaussc(l,k,j,-i)=akl
 925               endif
 926             enddo
 927           enddo
 928         enddo
 929         endif
 930       enddo
 931       close (irotam)
 932       if (lprint) then
 933         write (iout,'(/a)') 'Parameters of side-chain local geometry'
 934         do i=1,ntyp
 935           nlobi=nlob(i)
 936           if (nlobi.gt.0) then
 937           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
 938            ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
 939 !          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
 940 !          write (iout,'(a,f10.4,4(16x,f10.4))')
 941 !     &                             'Center  ',(bsc(j,i),j=1,nlobi)
 942 !          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
 943            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
 944                                    'log h',(bsc(j,i),j=1,nlobi)
 945            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
 946           'x',((censc(k,j,i),k=1,3),j=1,nlobi)
 947 !          write (iout,'(a)')
 948 !         do j=1,nlobi
 949 !           ind=0
 950 !           do k=1,3
 951 !             do l=1,k
 952 !              ind=ind+1
 953 !              blower(k,l,j)=gaussc(ind,j,i)
 954 !             enddo
 955 !           enddo
 956 !         enddo
 957           do k=1,3
 958             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
 959                        ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
 960           enddo
 961           endif
 962         enddo
 963       endif
 964 #else
 965 !
 966 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
 967 ! added by Urszula Kozlowska 07/11/2007
 968 !
 969       allocate(sc_parmin(65,ntyp))      !(maxsccoef,ntyp)
 970
 971       do i=1,ntyp
 972         read (irotam,*)
 973        if (i.eq.10) then
 974          read (irotam,*)
 975        else
 976          do j=1,65
 977            read(irotam,*) sc_parmin(j,i)
 978          enddo
 979        endif
 980       enddo
 981 #endif
 982       close(irotam)
 983 #ifdef CRYST_TOR
 984 !
 985 ! Read torsional parameters in old format
 986 !
 987       allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
 988
 989       read (itorp,*) ntortyp,nterm_old
 990       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
 991       read (itorp,*) (itortyp(i),i=1,ntyp)
 992
 993 !el from energy module--------
 994       allocate(v1(nterm_old,ntortyp,ntortyp))
 995       allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
 996 !el---------------------------
 997
 998       do i=1,ntortyp
 999         do j=1,ntortyp
1000           read (itorp,'(a)')
1001           do k=1,nterm_old
1002             read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
1003           enddo
1004         enddo
1005       enddo
1006       close (itorp)
1007       if (lprint) then
1008         write (iout,'(/a/)') 'Torsional constants:'
1009         do i=1,ntortyp
1010           do j=1,ntortyp
1011             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1012             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1013           enddo
1014         enddo
1015       endif
1016
1017
1018 #else
1019 !
1020 ! Read torsional parameters
1021 !
1022       allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1023
1024       read (itorp,*) ntortyp
1025       read (itorp,*) (itortyp(i),i=1,ntyp)
1026       write (iout,*) 'ntortyp',ntortyp
1027
1028 !el from energy module---------
1029       allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1030       allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1031
1032       allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1033       allocate(vlor2(maxlor,ntortyp,ntortyp))
1034       allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1035       allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1036
1037       allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1038       allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1039 !el---------------------------
1040       do iblock=1,2
1041         do i=-ntortyp,ntortyp
1042           do j=-ntortyp,ntortyp
1043             nterm(i,j,iblock)=0
1044             nlor(i,j,iblock)=0
1045           enddo
1046         enddo
1047       enddo
1048 !el---------------------------
1049
1050       do iblock=1,2
1051       do i=-ntyp,-1
1052        itortyp(i)=-itortyp(-i)
1053       enddo
1054 !      write (iout,*) 'ntortyp',ntortyp
1055       do i=0,ntortyp-1
1056         do j=-ntortyp+1,ntortyp-1
1057           read (itorp,*) nterm(i,j,iblock),&
1058                 nlor(i,j,iblock)
1059           nterm(-i,-j,iblock)=nterm(i,j,iblock)
1060           nlor(-i,-j,iblock)=nlor(i,j,iblock)
1061           v0ij=0.0d0
1062           si=-1.0d0
1063           do k=1,nterm(i,j,iblock)
1064             read (itorp,*) kk,v1(k,i,j,iblock),&
1065             v2(k,i,j,iblock)
1066             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1067             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1068             v0ij=v0ij+si*v1(k,i,j,iblock)
1069             si=-si
1070          enddo
1071           do k=1,nlor(i,j,iblock)
1072             read (itorp,*) kk,vlor1(k,i,j),&
1073               vlor2(k,i,j),vlor3(k,i,j)
1074             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1075           enddo
1076           v0(i,j,iblock)=v0ij
1077           v0(-i,-j,iblock)=v0ij
1078         enddo
1079       enddo
1080       enddo
1081       close (itorp)
1082       if (lprint) then
1083         do iblock=1,2 !el
1084         write (iout,'(/a/)') 'Torsional constants:'
1085         do i=1,ntortyp
1086           do j=1,ntortyp
1087             write (iout,*) 'ityp',i,' jtyp',j
1088             write (iout,*) 'Fourier constants'
1089             do k=1,nterm(i,j,iblock)
1090               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1091               v2(k,i,j,iblock)
1092             enddo
1093             write (iout,*) 'Lorenz constants'
1094             do k=1,nlor(i,j,iblock)
1095               write (iout,'(3(1pe15.5))') &
1096                vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1097             enddo
1098           enddo
1099         enddo
1100         enddo
1101       endif
1102 !
1103 ! 6/23/01 Read parameters for double torsionals
1104 !
1105 !el from energy module------------
1106       allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1107       allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1108 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1109       allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1110       allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1111         !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1112       allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1113       allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1114         !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1115 !---------------------------------
1116
1117       do iblock=1,2
1118       do i=0,ntortyp-1
1119         do j=-ntortyp+1,ntortyp-1
1120           do k=-ntortyp+1,ntortyp-1
1121             read (itordp,'(3a1)') t1,t2,t3
1122 !              write (iout,*) "OK onelett",
1123 !     &         i,j,k,t1,t2,t3
1124
1125             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1126               .or. t3.ne.toronelet(k)) then
1127               write (iout,*) "Error in double torsional parameter file",&
1128                i,j,k,t1,t2,t3
1129 #ifdef MPI
1130               call MPI_Finalize(Ierror)
1131 #endif
1132                stop "Error in double torsional parameter file"
1133             endif
1134           read (itordp,*) ntermd_1(i,j,k,iblock),&
1135                ntermd_2(i,j,k,iblock)
1136             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1137             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1138             read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,&
1139                ntermd_1(i,j,k,iblock))
1140             read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,&
1141                ntermd_1(i,j,k,iblock))
1142             read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,&
1143                ntermd_1(i,j,k,iblock))
1144             read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,&
1145                ntermd_1(i,j,k,iblock))
1146 ! Martix of D parameters for one dimesional foureir series
1147             do l=1,ntermd_1(i,j,k,iblock)
1148              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1149              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1150              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1151              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1152 !            write(iout,*) "whcodze" ,
1153 !     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1154             enddo
1155             read (itordp,*) ((v2c(l,m,i,j,k,iblock),&
1156                v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1157                v2s(m,l,i,j,k,iblock),&
1158                m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1159 ! Martix of D parameters for two dimesional fourier series
1160             do l=1,ntermd_2(i,j,k,iblock)
1161              do m=1,l-1
1162              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1163              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1164              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1165              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1166              enddo!m
1167             enddo!l
1168           enddo!k
1169         enddo!j
1170       enddo!i
1171       enddo!iblock
1172       if (lprint) then
1173       write (iout,*)
1174       write (iout,*) 'Constants for double torsionals'
1175       do iblock=1,2
1176       do i=0,ntortyp-1
1177         do j=-ntortyp+1,ntortyp-1
1178           do k=-ntortyp+1,ntortyp-1
1179             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1180               ' nsingle',ntermd_1(i,j,k,iblock),&
1181               ' ndouble',ntermd_2(i,j,k,iblock)
1182             write (iout,*)
1183             write (iout,*) 'Single angles:'
1184             do l=1,ntermd_1(i,j,k,iblock)
1185               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1186                  v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1187                  v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1188                  v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1189             enddo
1190             write (iout,*)
1191             write (iout,*) 'Pairs of angles:'
1192             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1193             do l=1,ntermd_2(i,j,k,iblock)
1194               write (iout,'(i5,20f10.5)') &
1195                l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1196             enddo
1197             write (iout,*)
1198            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1199             do l=1,ntermd_2(i,j,k,iblock)
1200               write (iout,'(i5,20f10.5)') &
1201                l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1202                (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1203             enddo
1204             write (iout,*)
1205           enddo
1206         enddo
1207       enddo
1208       enddo
1209       endif
1210 #endif
1211 !elwrite(iout,*) "parmread kontrol sc-bb"
1212 ! Read of Side-chain backbone correlation parameters
1213 ! Modified 11 May 2012 by Adasko
1214 !CC
1215 !
1216      read (isccor,*) nsccortyp
1217
1218      maxinter=3
1219 !c maxinter is maximum interaction sites
1220 !write(iout,*)"maxterm_sccor",maxterm_sccor
1221 !el from module energy-------------
1222       allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1223       allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1224       allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1225       allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp))   !(maxterm_sccor,20,20)
1226 !-----------------------------------
1227       allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1228 !-----------------------------------
1229       allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1230       allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1231                -nsccortyp:nsccortyp))
1232       allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1233                -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1234       allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1235                -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1236 !-----------------------------------
1237       do i=-nsccortyp,nsccortyp
1238         do j=-nsccortyp,nsccortyp
1239           nterm_sccor(j,i)=0
1240         enddo
1241       enddo
1242 !-----------------------------------
1243
1244       read (isccor,*) (isccortyp(i),i=1,ntyp)
1245       do i=-ntyp,-1
1246         isccortyp(i)=-isccortyp(-i)
1247       enddo
1248       iscprol=isccortyp(20)
1249 !      write (iout,*) 'ntortyp',ntortyp
1250 !      maxinter=3
1251 !c maxinter is maximum interaction sites
1252       do l=1,maxinter
1253       do i=1,nsccortyp
1254         do j=1,nsccortyp
1255           read (isccor,*) &
1256       nterm_sccor(i,j),nlor_sccor(i,j)
1257           v0ijsccor=0.0d0
1258           v0ijsccor1=0.0d0
1259           v0ijsccor2=0.0d0
1260           v0ijsccor3=0.0d0
1261           si=-1.0d0
1262           nterm_sccor(-i,j)=nterm_sccor(i,j)
1263           nterm_sccor(-i,-j)=nterm_sccor(i,j)
1264           nterm_sccor(i,-j)=nterm_sccor(i,j)
1265           do k=1,nterm_sccor(i,j)
1266             read (isccor,*) kk,v1sccor(k,l,i,j),&
1267             v2sccor(k,l,i,j)
1268             if (j.eq.iscprol) then
1269              if (i.eq.isccortyp(10)) then
1270              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1271              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1272              else
1273              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1274                               +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1275              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1276                               +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1277              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1278              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1279              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1280              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1281              endif
1282             else
1283              if (i.eq.isccortyp(10)) then
1284              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1285              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1286              else
1287                if (j.eq.isccortyp(10)) then
1288              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1289              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1290                else
1291              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1292              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1293              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1294              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1295              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1296              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1297                 endif
1298                endif
1299             endif
1300             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1301             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1302             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1303             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1304             si=-si
1305            enddo
1306           do k=1,nlor_sccor(i,j)
1307             read (isccor,*) kk,vlor1sccor(k,i,j),&
1308               vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1309             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1310       (1+vlor3sccor(k,i,j)**2)
1311           enddo
1312           v0sccor(l,i,j)=v0ijsccor
1313           v0sccor(l,-i,j)=v0ijsccor1
1314           v0sccor(l,i,-j)=v0ijsccor2
1315           v0sccor(l,-i,-j)=v0ijsccor3
1316           enddo
1317         enddo
1318       enddo
1319       close (isccor)
1320       if (lprint) then
1321         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
1322         do i=1,nsccortyp
1323           do j=1,nsccortyp
1324             write (iout,*) 'ityp',i,' jtyp',j
1325             write (iout,*) 'Fourier constants'
1326             do k=1,nterm_sccor(i,j)
1327               write (iout,'(2(1pe15.5))') &
1328          (v1sccor(k,l,i,j),v2sccor(k,l,i,j),l=1,maxinter)
1329             enddo
1330             write (iout,*) 'Lorenz constants'
1331             do k=1,nlor_sccor(i,j)
1332               write (iout,'(3(1pe15.5))') &
1333                vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1334             enddo
1335           enddo
1336         enddo
1337       endif
1338 !
1339 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1340 !         interaction energy of the Gly, Ala, and Pro prototypes.
1341 !
1342       read (ifourier,*) nloctyp
1343 !el write(iout,*)"nloctyp",nloctyp
1344 !el from module energy-------
1345       allocate(b1(2,-nloctyp-1:nloctyp+1))      !(2,-maxtor:maxtor)
1346       allocate(b2(2,-nloctyp-1:nloctyp+1))      !(2,-maxtor:maxtor)
1347       allocate(b1tilde(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1348       allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1349       allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1350       allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1351       allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1352       allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1353       do i=1,2
1354         do ii=-nloctyp-1,nloctyp+1
1355           b1(i,ii)=0.0d0
1356           b2(i,ii)=0.0d0
1357           b1tilde(i,ii)=0.0d0
1358           do j=1,2
1359             cc(j,i,ii)=0.0d0
1360             dd(j,i,ii)=0.0d0
1361             ee(j,i,ii)=0.0d0
1362             ctilde(j,i,ii)=0.0d0
1363             dtilde(j,i,ii)=0.0d0
1364           enddo
1365         enddo
1366       enddo
1367 !--------------------------------
1368       allocate(b(13,0:nloctyp))
1369
1370       do i=0,nloctyp-1
1371         read (ifourier,*)
1372         read (ifourier,*) (b(ii,i),ii=1,13)
1373         if (lprint) then
1374         write (iout,*) 'Type',i
1375         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1376         endif
1377         B1(1,i)  = b(3,i)
1378         B1(2,i)  = b(5,i)
1379         B1(1,-i) = b(3,i)
1380         B1(2,-i) = -b(5,i)
1381 !        b1(1,i)=0.0d0
1382 !        b1(2,i)=0.0d0
1383         B1tilde(1,i) = b(3,i)
1384         B1tilde(2,i) =-b(5,i)
1385         B1tilde(1,-i) =-b(3,i)
1386         B1tilde(2,-i) =b(5,i)
1387 !        b1tilde(1,i)=0.0d0
1388 !        b1tilde(2,i)=0.0d0
1389         B2(1,i)  = b(2,i)
1390         B2(2,i)  = b(4,i)
1391         B2(1,-i)  =b(2,i)
1392         B2(2,-i)  =-b(4,i)
1393
1394 !        b2(1,i)=0.0d0
1395 !        b2(2,i)=0.0d0
1396         CC(1,1,i)= b(7,i)
1397         CC(2,2,i)=-b(7,i)
1398         CC(2,1,i)= b(9,i)
1399         CC(1,2,i)= b(9,i)
1400         CC(1,1,-i)= b(7,i)
1401         CC(2,2,-i)=-b(7,i)
1402         CC(2,1,-i)=-b(9,i)
1403         CC(1,2,-i)=-b(9,i)
1404 !        CC(1,1,i)=0.0d0
1405 !        CC(2,2,i)=0.0d0
1406 !        CC(2,1,i)=0.0d0
1407 !        CC(1,2,i)=0.0d0
1408         Ctilde(1,1,i)=b(7,i)
1409         Ctilde(1,2,i)=b(9,i)
1410         Ctilde(2,1,i)=-b(9,i)
1411         Ctilde(2,2,i)=b(7,i)
1412         Ctilde(1,1,-i)=b(7,i)
1413         Ctilde(1,2,-i)=-b(9,i)
1414         Ctilde(2,1,-i)=b(9,i)
1415         Ctilde(2,2,-i)=b(7,i)
1416
1417 !        Ctilde(1,1,i)=0.0d0
1418 !        Ctilde(1,2,i)=0.0d0
1419 !        Ctilde(2,1,i)=0.0d0
1420 !        Ctilde(2,2,i)=0.0d0
1421         DD(1,1,i)= b(6,i)
1422         DD(2,2,i)=-b(6,i)
1423         DD(2,1,i)= b(8,i)
1424         DD(1,2,i)= b(8,i)
1425         DD(1,1,-i)= b(6,i)
1426         DD(2,2,-i)=-b(6,i)
1427         DD(2,1,-i)=-b(8,i)
1428         DD(1,2,-i)=-b(8,i)
1429 !        DD(1,1,i)=0.0d0
1430 !        DD(2,2,i)=0.0d0
1431 !        DD(2,1,i)=0.0d0
1432 !        DD(1,2,i)=0.0d0
1433         Dtilde(1,1,i)=b(6,i)
1434         Dtilde(1,2,i)=b(8,i)
1435         Dtilde(2,1,i)=-b(8,i)
1436         Dtilde(2,2,i)=b(6,i)
1437         Dtilde(1,1,-i)=b(6,i)
1438         Dtilde(1,2,-i)=-b(8,i)
1439         Dtilde(2,1,-i)=b(8,i)
1440         Dtilde(2,2,-i)=b(6,i)
1441
1442 !        Dtilde(1,1,i)=0.0d0
1443 !        Dtilde(1,2,i)=0.0d0
1444 !        Dtilde(2,1,i)=0.0d0
1445 !        Dtilde(2,2,i)=0.0d0
1446         EE(1,1,i)= b(10,i)+b(11,i)
1447         EE(2,2,i)=-b(10,i)+b(11,i)
1448         EE(2,1,i)= b(12,i)-b(13,i)
1449         EE(1,2,i)= b(12,i)+b(13,i)
1450         EE(1,1,-i)= b(10,i)+b(11,i)
1451         EE(2,2,-i)=-b(10,i)+b(11,i)
1452         EE(2,1,-i)=-b(12,i)+b(13,i)
1453         EE(1,2,-i)=-b(12,i)-b(13,i)
1454
1455 !        ee(1,1,i)=1.0d0
1456 !        ee(2,2,i)=1.0d0
1457 !        ee(2,1,i)=0.0d0
1458 !        ee(1,2,i)=0.0d0
1459 !        ee(2,1,i)=ee(1,2,i)
1460
1461       enddo
1462       if (lprint) then
1463       do i=1,nloctyp
1464         write (iout,*) 'Type',i
1465         write (iout,*) 'B1'
1466 !        write (iout,'(f10.5)') B1(:,i)
1467         write(iout,*) B1(1,i),B1(2,i)
1468         write (iout,*) 'B2'
1469 !        write (iout,'(f10.5)') B2(:,i)
1470         write(iout,*) B2(1,i),B2(2,i)
1471         write (iout,*) 'CC'
1472         do j=1,2
1473           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1474         enddo
1475         write(iout,*) 'DD'
1476         do j=1,2
1477           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1478         enddo
1479         write(iout,*) 'EE'
1480         do j=1,2
1481           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1482         enddo
1483       enddo
1484       endif
1485 !
1486 ! Read electrostatic-interaction parameters
1487 !
1488       if (lprint) then
1489         write (iout,'(/a)') 'Electrostatic interaction constants:'
1490         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
1491                   'IT','JT','APP','BPP','AEL6','AEL3'
1492       endif
1493       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
1494       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
1495       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
1496       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
1497       close (ielep)
1498       do i=1,2
1499         do j=1,2
1500         rri=rpp(i,j)**6
1501         app (i,j)=epp(i,j)*rri*rri
1502         bpp (i,j)=-2.0D0*epp(i,j)*rri
1503         ael6(i,j)=elpp6(i,j)*4.2D0**6
1504         ael3(i,j)=elpp3(i,j)*4.2D0**3
1505         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
1506                           ael6(i,j),ael3(i,j)
1507         enddo
1508       enddo
1509 !
1510 ! Read side-chain interaction parameters.
1511 !
1512 !el from module energy - COMMON.INTERACT-------
1513       allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
1514       allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
1515       allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
1516       allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
1517       allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
1518       do i=1,ntyp
1519         do j=1,ntyp
1520           augm(i,j)=0.0D0
1521         enddo
1522         chip(i)=0.0D0
1523         alp(i)=0.0D0
1524         sigma0(i)=0.0D0
1525         sigii(i)=0.0D0
1526         rr0(i)=0.0D0
1527       enddo
1528 !--------------------------------
1529
1530       read (isidep,*) ipot,expon
1531 !el      if (ipot.lt.1 .or. ipot.gt.5) then
1532 !        write (iout,'(2a)') 'Error while reading SC interaction',&
1533 !                     'potential file - unknown potential type.'
1534 !        stop
1535 !wl      endif
1536       expon2=expon/2
1537       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
1538        ', exponents are ',expon,2*expon
1539 !      goto (10,20,30,30,40) ipot
1540       select case(ipot)
1541 !----------------------- LJ potential ---------------------------------
1542        case (1)
1543 !   10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1544         read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
1545         if (lprint) then
1546           write (iout,'(/a/)') 'Parameters of the LJ potential:'
1547           write (iout,'(a/)') 'The epsilon array:'
1548           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1549           write (iout,'(/a)') 'One-body parameters:'
1550           write (iout,'(a,4x,a)') 'residue','sigma'
1551           write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1552         endif
1553 !      goto 50
1554 !----------------------- LJK potential --------------------------------
1555        case (2)
1556 !   20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1557         read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1558           (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1559         if (lprint) then
1560           write (iout,'(/a/)') 'Parameters of the LJK potential:'
1561           write (iout,'(a/)') 'The epsilon array:'
1562           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1563           write (iout,'(/a)') 'One-body parameters:'
1564           write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1565           write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
1566                 i=1,ntyp)
1567         endif
1568 !      goto 50
1569 !---------------------- GB or BP potential -----------------------------
1570        case (3:4)
1571 !   30 do i=1,ntyp
1572         do i=1,ntyp
1573          read (isidep,*)(eps(i,j),j=i,ntyp)
1574         enddo
1575         read (isidep,*)(sigma0(i),i=1,ntyp)
1576         read (isidep,*)(sigii(i),i=1,ntyp)
1577         read (isidep,*)(chip(i),i=1,ntyp)
1578         read (isidep,*)(alp(i),i=1,ntyp)
1579 ! For the GB potential convert sigma'**2 into chi'
1580         if (ipot.eq.4) then
1581           do i=1,ntyp
1582             chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1583           enddo
1584         endif
1585         if (lprint) then
1586           write (iout,'(/a/)') 'Parameters of the BP potential:'
1587           write (iout,'(a/)') 'The epsilon array:'
1588           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1589           write (iout,'(/a)') 'One-body parameters:'
1590           write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',&
1591                '    chip  ','    alph  '
1592           write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
1593                            chip(i),alp(i),i=1,ntyp)
1594         endif
1595 !      goto 50
1596 !--------------------- GBV potential -----------------------------------
1597        case (5)
1598 !   40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1599         read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1600           (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
1601         (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
1602         if (lprint) then
1603           write (iout,'(/a/)') 'Parameters of the GBV potential:'
1604           write (iout,'(a/)') 'The epsilon array:'
1605           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1606           write (iout,'(/a)') 'One-body parameters:'
1607           write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',&
1608             's||/s_|_^2','    chip  ','    alph  '
1609           write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
1610                  sigii(i),chip(i),alp(i),i=1,ntyp)
1611         endif
1612        case default
1613         write (iout,'(2a)') 'Error while reading SC interaction',&
1614                      'potential file - unknown potential type.'
1615         stop
1616 !   50 continue
1617       end select
1618 !      continue
1619       close (isidep)
1620 !-----------------------------------------------------------------------
1621 ! Calculate the "working" parameters of SC interactions.
1622
1623 !el from module energy - COMMON.INTERACT-------
1624       allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
1625       allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
1626       do i=1,ntyp1
1627         do j=1,ntyp1
1628           aa(i,j)=0.0D0
1629           bb(i,j)=0.0D0
1630           chi(i,j)=0.0D0
1631           sigma(i,j)=0.0D0
1632           r0(i,j)=0.0D0
1633         enddo
1634       enddo
1635 !--------------------------------
1636
1637       do i=2,ntyp
1638         do j=1,i-1
1639           eps(i,j)=eps(j,i)
1640         enddo
1641       enddo
1642       do i=1,ntyp
1643         do j=i,ntyp
1644           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1645           sigma(j,i)=sigma(i,j)
1646           rs0(i,j)=dwa16*sigma(i,j)
1647           rs0(j,i)=rs0(i,j)
1648         enddo
1649       enddo
1650       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
1651        'Working parameters of the SC interactions:',&
1652        '     a    ','     b    ','   augm   ','  sigma ','   r0   ',&
1653        '  chi1   ','   chi2   '
1654       do i=1,ntyp
1655         do j=i,ntyp
1656           epsij=eps(i,j)
1657           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1658             rrij=sigma(i,j)
1659           else
1660             rrij=rr0(i)+rr0(j)
1661           endif
1662           r0(i,j)=rrij
1663           r0(j,i)=rrij
1664           rrij=rrij**expon
1665           epsij=eps(i,j)
1666           sigeps=dsign(1.0D0,epsij)
1667           epsij=dabs(epsij)
1668           aa(i,j)=epsij*rrij*rrij
1669           bb(i,j)=-sigeps*epsij*rrij
1670           aa(j,i)=aa(i,j)
1671           bb(j,i)=bb(i,j)
1672           if (ipot.gt.2) then
1673             sigt1sq=sigma0(i)**2
1674             sigt2sq=sigma0(j)**2
1675             sigii1=sigii(i)
1676             sigii2=sigii(j)
1677             ratsig1=sigt2sq/sigt1sq
1678             ratsig2=1.0D0/ratsig1
1679             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1680             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1681             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1682           else
1683             rsum_max=sigma(i,j)
1684           endif
1685 !         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1686             sigmaii(i,j)=rsum_max
1687             sigmaii(j,i)=rsum_max
1688 !         else
1689 !           sigmaii(i,j)=r0(i,j)
1690 !           sigmaii(j,i)=r0(i,j)
1691 !         endif
1692 !d        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
1693           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
1694             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1695             augm(i,j)=epsij*r_augm**(2*expon)
1696 !           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1697             augm(j,i)=augm(i,j)
1698           else
1699             augm(i,j)=0.0D0
1700             augm(j,i)=0.0D0
1701           endif
1702           if (lprint) then
1703             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')  &
1704             restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),&
1705             sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1706           endif
1707         enddo
1708       enddo
1709
1710       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
1711       do i=1,ntyp
1712         do j=1,2
1713           bad(i,j)=0.0D0
1714         enddo
1715       enddo
1716 #ifdef CLUSTER
1717 !
1718 ! Define the SC-p interaction constants
1719 !
1720       do i=1,20
1721         do j=1,2
1722           eps_scp(i,j)=-1.5d0
1723           rscp(i,j)=4.0d0
1724         enddo
1725       enddo
1726 #endif
1727
1728 !elwrite(iout,*) "parmread kontrol before oldscp"
1729 !
1730 ! Define the SC-p interaction constants
1731 !
1732 #ifdef OLDSCP
1733       do i=1,20
1734 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
1735 ! helix formation)
1736 !       aad(i,1)=0.3D0*4.0D0**12
1737 ! Following line for constants currently implemented
1738 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
1739         aad(i,1)=1.5D0*4.0D0**12
1740 !       aad(i,1)=0.17D0*5.6D0**12
1741         aad(i,2)=aad(i,1)
1742 ! "Soft" SC-p repulsion
1743         bad(i,1)=0.0D0
1744 ! Following line for constants currently implemented
1745 !       aad(i,1)=0.3D0*4.0D0**6
1746 ! "Hard" SC-p repulsion
1747         bad(i,1)=3.0D0*4.0D0**6
1748 !       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
1749         bad(i,2)=bad(i,1)
1750 !       aad(i,1)=0.0D0
1751 !       aad(i,2)=0.0D0
1752 !       bad(i,1)=1228.8D0
1753 !       bad(i,2)=1228.8D0
1754       enddo
1755 #else
1756 !
1757 ! 8/9/01 Read the SC-p interaction constants from file
1758 !
1759       do i=1,ntyp
1760         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
1761       enddo
1762       do i=1,ntyp
1763         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1764         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1765         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1766         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1767       enddo
1768
1769       if (lprint) then
1770         write (iout,*) "Parameters of SC-p interactions:"
1771         do i=1,20
1772           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
1773            eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1774         enddo
1775       endif
1776 #endif
1777 !
1778 ! Define the constants of the disulfide bridge
1779 !
1780       ebr=-5.50D0
1781 !
1782 ! Old arbitrary potential - commented out.
1783 !
1784 !      dbr= 4.20D0
1785 !      fbr= 3.30D0
1786 !
1787 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1788 ! energy surface of diethyl disulfide.
1789 ! A. Liwo and U. Kozlowska, 11/24/03
1790 !
1791       D0CM = 3.78d0
1792       AKCM = 15.1d0
1793       AKTH = 11.0d0
1794       AKCT = 12.0d0
1795       V1SS =-1.08d0
1796       V2SS = 7.61d0
1797       V3SS = 13.7d0
1798
1799       if (lprint) then
1800       write (iout,'(/a)') "Disulfide bridge parameters:"
1801       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1802       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1803       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1804       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
1805        ' v3ss:',v3ss
1806       endif
1807       return
1808       end subroutine parmread
1809 #ifndef CLUSTER
1810 !-----------------------------------------------------------------------------
1811 ! mygetenv.F
1812 !-----------------------------------------------------------------------------
1813       subroutine mygetenv(string,var)
1814 !
1815 ! Version 1.0
1816 !
1817 ! This subroutine passes the environmental variables to FORTRAN program.
1818 ! If the flags -DMYGETENV and -DMPI are not for compilation, it calls the
1819 ! standard FORTRAN GETENV subroutine. If both flags are set, the subroutine
1820 ! reads the environmental variables from $HOME/.env
1821 !
1822 ! Usage: As for the standard FORTRAN GETENV subroutine.
1823 !
1824 ! Purpose: some versions/installations of MPI do not transfer the environmental
1825 ! variables to slave processors, if these variables are set in the shell script
1826 ! from which mpirun is called.
1827 !
1828 ! A.Liwo, 7/29/01
1829 !
1830 #ifdef MPI
1831       use MPI_data
1832       include "mpif.h"
1833 #endif
1834 !      implicit none
1835       character*(*) :: string,var
1836 #if defined(MYGETENV) && defined(MPI)
1837 !      include "DIMENSIONS.ZSCOPT"
1838 !      include "mpif.h"
1839 !      include "COMMON.MPI"
1840 !el      character*360 ucase
1841 !el      external ucase
1842       character(len=360) :: string1(360),karta
1843       character(len=240) :: home
1844       integer i,n !,ilen
1845 !el      external ilen
1846       call getenv("HOME",home)
1847       open(99,file=home(:ilen(home))//"/.env",status="OLD",err=112)
1848       do while (.true.)
1849         read (99,end=111,err=111,'(a)') karta
1850         do i=1,80
1851           string1(i)=" "
1852         enddo
1853         call split_string(karta,string1,80,n)
1854         if (ucase(string1(1)(:ilen(string1(1)))).eq."SETENV" .and. &
1855          string1(2)(:ilen(string1(2))).eq.string(:ilen(string)) ) then
1856            var=string1(3)
1857            print *,"Processor",me,": ",var(:ilen(var)),&
1858             " assigned to ",string(:ilen(string))
1859            close(99)
1860            return
1861         endif
1862       enddo
1863  111  print *,"Environment variable ",string(:ilen(string))," not set."
1864       close(99)
1865       return
1866  112  print *,"Error opening environment file!"
1867 #else
1868       call getenv(string,var)
1869 #endif
1870       return
1871       end subroutine mygetenv
1872 !-----------------------------------------------------------------------------
1873 ! readrtns.F
1874 !-----------------------------------------------------------------------------
1875       subroutine read_general_data(*)
1876
1877       use control_data, only:indpdb,symetr
1878       use energy_data, only:distchainmax
1879 !      implicit none
1880 !      include "DIMENSIONS"
1881 !      include "DIMENSIONS.ZSCOPT"
1882 !      include "DIMENSIONS.FREE"
1883 !      include "COMMON.TORSION"
1884 !      include "COMMON.INTERACT"
1885 !      include "COMMON.IOUNITS"
1886 !      include "COMMON.TIME1"
1887 !      include "COMMON.PROT"
1888 !      include "COMMON.PROTFILES"
1889 !      include "COMMON.CHAIN"
1890 !      include "COMMON.NAMES"
1891 !      include "COMMON.FFIELD"
1892 !      include "COMMON.ENEPS"
1893 !      include "COMMON.WEIGHTS"
1894 !      include "COMMON.FREE"
1895 !      include "COMMON.CONTROL"
1896 !      include "COMMON.ENERGIES"
1897       character(len=800) :: controlcard
1898       integer :: i,j,k,ii,n_ene_found
1899       integer :: ind,itype1,itype2,itypf,itypsc,itypp
1900 !el      integer ilen
1901 !el      external ilen
1902 !el      character*16 ucase
1903       character(len=16) :: key
1904 !el      external ucase
1905       call card_concat(controlcard,.true.)
1906       call readi(controlcard,"N_ENE",n_eneW,max_eneW)
1907       if (n_eneW.gt.max_eneW) then
1908         write (iout,*) "Error: parameter out of range: N_ENE",n_eneW,&
1909           max_eneW
1910         return 1
1911       endif
1912       call readi(controlcard,"NPARMSET",nparmset,1)
1913 !elwrite(iout,*)"in read_gen data"
1914       separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
1915       call readi(controlcard,"IPARMPRINT",iparmprint,1)
1916       write (iout,*) "PARMPRINT",iparmprint
1917       if (nparmset.gt.max_parm) then
1918         write (iout,*) "Error: parameter out of range: NPARMSET",&
1919           nparmset, Max_Parm
1920         return 1
1921       endif
1922 !elwrite(iout,*)"in read_gen data"
1923       call readi(controlcard,"MAXIT",maxit,5000)
1924       call reada(controlcard,"FIMIN",fimin,1.0d-3)
1925       call readi(controlcard,"ENSEMBLES",ensembles,0)
1926       hamil_rep=index(controlcard,"HAMIL_REP").gt.0
1927       write (iout,*) "Number of energy parameter sets",nparmset
1928       allocate(isampl(nparmset))
1929       call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
1930       write (iout,*) "MaxSlice",MaxSlice
1931       call readi(controlcard,"NSLICE",nslice,1)
1932 !elwrite(iout,*)"in read_gen data"
1933       call flush(iout)
1934       if (nslice.gt.MaxSlice) then
1935         write (iout,*) "Error: parameter out of range: NSLICE",nslice,&
1936           MaxSlice
1937         return 1
1938       endif
1939       write (iout,*) "Frequency of storing conformations",&
1940        (isampl(i),i=1,nparmset)
1941       write (iout,*) "Maxit",maxit," Fimin",fimin
1942       call readi(controlcard,"NQ",nQ,1)
1943       if (nQ.gt.MaxQ) then
1944         write (iout,*) "Error: parameter out of range: NQ",nq,&
1945           maxq
1946         return 1
1947       endif
1948       indpdb=0
1949       if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
1950       call reada(controlcard,"DELTA",delta,1.0d-2)
1951       call readi(controlcard,"EINICHECK",einicheck,2)
1952       call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
1953       call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
1954       call readi(controlcard,"RESCALE",rescale_modeW,1)
1955       check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
1956       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
1957       call readi(controlcard,'SYM',symetr,1)
1958       write (iout,*) "DISTCHAINMAX",distchainmax
1959       write (iout,*) "delta",delta
1960       write (iout,*) "einicheck",einicheck
1961       write (iout,*) "rescale_mode",rescale_modeW
1962       call flush(iout)
1963       bxfile=index(controlcard,"BXFILE").gt.0
1964       cxfile=index(controlcard,"CXFILE").gt.0
1965       if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile) &
1966        bxfile=.true.
1967       histfile=index(controlcard,"HISTFILE").gt.0
1968       histout=index(controlcard,"HISTOUT").gt.0
1969       entfile=index(controlcard,"ENTFILE").gt.0
1970       zscfile=index(controlcard,"ZSCFILE").gt.0
1971       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1972       write (iout,*) "with_dihed_constr ",with_dihed_constr
1973       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
1974       return
1975       end subroutine read_general_data
1976 !------------------------------------------------------------------------------
1977       subroutine read_efree(*)
1978 !
1979 ! Read molecular data
1980 !
1981 !      implicit none
1982 !      include 'DIMENSIONS'
1983 !      include 'DIMENSIONS.ZSCOPT'
1984 !      include 'DIMENSIONS.COMPAR'
1985 !      include 'DIMENSIONS.FREE'
1986 !      include 'COMMON.IOUNITS'
1987 !      include 'COMMON.TIME1'
1988 !      include 'COMMON.SBRIDGE'
1989 !      include 'COMMON.CONTROL'
1990 !      include 'COMMON.CHAIN'
1991 !      include 'COMMON.HEADER'
1992 !      include 'COMMON.GEO'
1993 !      include 'COMMON.FREE'
1994       character(len=320) :: controlcard !,ucase
1995       integer :: iparm,ib,i,j,npars
1996 !el      integer ilen
1997 !el      external ilen
1998
1999       if (hamil_rep) then
2000         npars=1
2001       else
2002         npars=nParmSet
2003       endif
2004
2005 !      call alloc_wham_arrays
2006 !      allocate(nT_h(nParmSet))
2007 !      allocate(replica(nParmSet))
2008 !      allocate(umbrella(nParmSet))
2009 !      allocate(read_iset(nParmSet))
2010 !      allocate(nT_h(nParmSet))
2011
2012       do iparm=1,npars
2013
2014       call card_concat(controlcard,.true.)
2015       call readi(controlcard,'NT',nT_h(iparm),1)
2016       write (iout,*) "iparm",iparm," nt",nT_h(iparm)
2017       call flush(iout)
2018       if (nT_h(iparm).gt.MaxT_h) then
2019         write (iout,*)  "Error: parameter out of range: NT",nT_h(iparm),&
2020           MaxT_h
2021         return 1
2022       endif
2023       replica(iparm)=index(controlcard,"REPLICA").gt.0
2024       umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
2025       read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
2026       write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",&
2027         replica(iparm)," umbrella ",umbrella(iparm),&
2028         " read_iset",read_iset(iparm)
2029       call flush(iout)
2030       do ib=1,nT_h(iparm)
2031         call card_concat(controlcard,.true.)
2032         call readi(controlcard,'NR',nR(ib,iparm),1)
2033         if (umbrella(iparm)) then
2034           nRR(ib,iparm)=1
2035         else
2036           nRR(ib,iparm)=nR(ib,iparm)
2037         endif
2038         if (nR(ib,iparm).gt.MaxR) then
2039           write (iout,*)  "Error: parameter out of range: NR",&
2040             nR(ib,iparm),MaxR
2041         return 1
2042         endif
2043         call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
2044         beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
2045         call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),&
2046           0.0d0)
2047         do i=1,nR(ib,iparm)
2048           call card_concat(controlcard,.true.)
2049           call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,&
2050             100.0d0)
2051           call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,&
2052             0.0d0)
2053         enddo
2054       enddo
2055       do ib=1,nT_h(iparm)
2056         write (iout,*) "ib",ib," beta_h",&
2057           1.0d0/(0.001987*beta_h(ib,iparm))
2058         write (iout,*) "nR",nR(ib,iparm)
2059         write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
2060         do i=1,nR(ib,iparm)
2061           write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),&
2062             "q0",(q0(j,i,ib,iparm),j=1,nQ)
2063         enddo
2064         call flush(iout)
2065       enddo
2066
2067       enddo
2068
2069       if (hamil_rep) then
2070
2071        do iparm=2,nParmSet
2072           nT_h(iparm)=nT_h(1)
2073          do ib=1,nT_h(iparm)
2074            nR(ib,iparm)=nR(ib,1)
2075            if (umbrella(iparm)) then
2076              nRR(ib,iparm)=1
2077            else
2078              nRR(ib,iparm)=nR(ib,1)
2079            endif
2080            beta_h(ib,iparm)=beta_h(ib,1)
2081            do i=1,nR(ib,iparm)
2082              f(i,ib,iparm)=f(i,ib,1)
2083              do j=1,nQ
2084                KH(j,i,ib,iparm)=KH(j,i,ib,1)
2085                Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
2086              enddo
2087            enddo
2088            replica(iparm)=replica(1)
2089            umbrella(iparm)=umbrella(1)
2090            read_iset(iparm)=read_iset(1)
2091          enddo
2092        enddo
2093
2094       endif
2095
2096       return
2097       end subroutine read_efree
2098 !-----------------------------------------------------------------------------
2099       subroutine read_protein_data(*)
2100 !      implicit none
2101 !      include "DIMENSIONS"
2102 !      include "DIMENSIONS.ZSCOPT"
2103 !      include "DIMENSIONS.FREE"
2104 #ifdef MPI
2105       use MPI_data
2106       include "mpif.h"
2107       integer :: IERROR,ERRCODE!,STATUS(MPI_STATUS_SIZE)
2108 !      include "COMMON.MPI"
2109 #endif
2110 !      include "COMMON.CHAIN"
2111 !      include "COMMON.IOUNITS"
2112 !      include "COMMON.PROT"
2113 !      include "COMMON.PROTFILES"
2114 !      include "COMMON.NAMES"
2115 !      include "COMMON.FREE"
2116 !      include "COMMON.OBCINKA"
2117       character(len=64) :: nazwa
2118       character(len=16000) :: controlcard
2119       integer :: i,ii,ib,iR,iparm,nthr,npars !,ilen,iroof
2120 !el      external ilen,iroof
2121       if (hamil_rep) then
2122         npars=1
2123       else
2124         npars=nparmset
2125       endif
2126
2127       do iparm=1,npars
2128
2129 ! Read names of files with conformation data.
2130       if (replica(iparm)) then
2131         nthr = 1
2132       else
2133         nthr = nT_h(iparm)
2134       endif
2135       do ib=1,nthr
2136       do ii=1,nRR(ib,iparm)
2137       write (iout,*) "Parameter set",iparm," temperature",ib,&
2138        " window",ii
2139       call flush(iout)
2140       call card_concat(controlcard,.true.)
2141       write (iout,*) controlcard(:ilen(controlcard))
2142       call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
2143       call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
2144       call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
2145       call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
2146       call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),&
2147        maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
2148       call reada(controlcard,"TIME_START",&
2149         time_start_collect(ii,ib,iparm),0.0d0)
2150       call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),&
2151         1.0d10)
2152       write (iout,*) "rec_start",rec_start(ii,ib,iparm),&
2153        " rec_end",rec_end(ii,ib,iparm)
2154       write (iout,*) "time_start",time_start_collect(ii,ib,iparm),&
2155        " time_end",time_end_collect(ii,ib,iparm)
2156       call flush(iout)
2157       if (replica(iparm)) then
2158         call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
2159         write (iout,*) "Number of trajectories",totraj(ii,iparm)
2160         call flush(iout)
2161       endif
2162       if (nfile_bin(ii,ib,iparm).lt.2 &
2163           .and. nfile_asc(ii,ib,iparm).eq.0 &
2164           .and. nfile_cx(ii,ib,iparm).eq.0) then
2165         write (iout,*) "Error - no action specified!"
2166         return 1
2167       endif
2168       if (nfile_bin(ii,ib,iparm).gt.0) then
2169         call card_concat(controlcard,.false.)
2170         call split_string(controlcard,protfiles(1,1,ii,ib,iparm),&
2171          maxfile_prot,nfile_bin(ii,ib,iparm))
2172 #ifdef DEBUG
2173         write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
2174         write(iout,*) (protfiles(i,1,ii,ib,iparm),&
2175           i=1,nfile_bin(ii,ib,iparm))
2176 #endif
2177       endif
2178       if (nfile_asc(ii,ib,iparm).gt.0) then
2179         call card_concat(controlcard,.false.)
2180         call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
2181          maxfile_prot,nfile_asc(ii,ib,iparm))
2182 #ifdef DEBUG
2183         write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
2184         write(iout,*) (protfiles(i,2,ii,ib,iparm),&
2185           i=1,nfile_asc(ii,ib,iparm))
2186 #endif
2187       else if (nfile_cx(ii,ib,iparm).gt.0) then
2188         call card_concat(controlcard,.false.)
2189         call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
2190          maxfile_prot,nfile_cx(ii,ib,iparm))
2191 #ifdef DEBUG
2192         write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
2193         write(iout,*) (protfiles(i,2,ii,ib,iparm),&
2194          i=1,nfile_cx(ii,ib,iparm))
2195 #endif
2196       endif
2197       call flush(iout)
2198       enddo
2199       enddo
2200
2201       enddo
2202       return
2203       end subroutine read_protein_data
2204 !-------------------------------------------------------------------------------
2205       subroutine readsss(rekord,lancuch,wartosc,default)
2206 !      implicit none
2207       character*(*) :: rekord,lancuch,wartosc,default
2208       character(len=80) :: aux
2209       integer :: lenlan,lenrec,iread,ireade
2210 !el      external ilen
2211 !el      logical iblnk
2212 !el      external iblnk
2213       lenlan=ilen(lancuch)
2214       lenrec=ilen(rekord)
2215       iread=index(rekord,lancuch(:lenlan)//"=")
2216 !      print *,"rekord",rekord," lancuch",lancuch
2217 !      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
2218       if (iread.eq.0) then
2219         wartosc=default
2220         return
2221       endif
2222       iread=iread+lenlan+1
2223 !      print *,"iread",iread
2224 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2225       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
2226         iread=iread+1
2227 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2228       enddo
2229 !      print *,"iread",iread
2230       if (iread.gt.lenrec) then
2231          wartosc=default
2232         return
2233       endif
2234       ireade=iread+1
2235 !      print *,"ireade",ireade
2236       do while (ireade.lt.lenrec .and. &
2237          .not.iblnk(rekord(ireade:ireade)))
2238         ireade=ireade+1
2239       enddo
2240       wartosc=rekord(iread:ireade)
2241       return
2242       end subroutine readsss
2243 !----------------------------------------------------------------------------
2244       subroutine multreads(rekord,lancuch,tablica,dim,default)
2245 !      implicit none
2246       integer :: dim,i
2247       character*(*) rekord,lancuch,tablica(dim),default
2248       character(len=80) :: aux
2249       integer :: lenlan,lenrec,iread,ireade
2250 !el      external ilen
2251 !el      logical iblnk
2252 !el      external iblnk
2253       do i=1,dim
2254         tablica(i)=default
2255       enddo
2256       lenlan=ilen(lancuch)
2257       lenrec=ilen(rekord)
2258       iread=index(rekord,lancuch(:lenlan)//"=")
2259 !      print *,"rekord",rekord," lancuch",lancuch
2260 !      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
2261       if (iread.eq.0) return
2262       iread=iread+lenlan+1
2263       do i=1,dim
2264 !      print *,"iread",iread
2265 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2266       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
2267         iread=iread+1
2268 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2269       enddo
2270 !      print *,"iread",iread
2271       if (iread.gt.lenrec) return
2272       ireade=iread+1
2273 !      print *,"ireade",ireade
2274       do while (ireade.lt.lenrec .and. &
2275          .not.iblnk(rekord(ireade:ireade)))
2276         ireade=ireade+1
2277       enddo
2278       tablica(i)=rekord(iread:ireade)
2279       iread=ireade+1
2280       enddo
2281       end subroutine multreads
2282 !----------------------------------------------------------------------------
2283       subroutine split_string(rekord,tablica,dim,nsub)
2284 !      implicit none
2285       integer :: dim,nsub,i,ii,ll,kk
2286       character*(*) tablica(dim)
2287       character*(*) rekord
2288 !el      integer ilen
2289 !el      external ilen
2290       do i=1,dim
2291         tablica(i)=" "
2292       enddo
2293       ii=1
2294       ll = ilen(rekord)
2295       nsub=0
2296       do i=1,dim
2297 ! Find the start of term name
2298         kk = 0
2299         do while (ii.le.ll .and. rekord(ii:ii).eq." ")
2300           ii = ii+1
2301         enddo
2302 ! Parse the name into TABLICA(i) until blank found
2303         do while (ii.le.ll .and. rekord(ii:ii).ne." ")
2304           kk = kk+1
2305           tablica(i)(kk:kk)=rekord(ii:ii)
2306           ii = ii+1
2307         enddo
2308         if (kk.gt.0) nsub=nsub+1
2309         if (ii.gt.ll) return
2310       enddo
2311       return
2312       end subroutine split_string
2313 !--------------------------------------------------------------------------------
2314 ! readrtns_compar.F
2315 !--------------------------------------------------------------------------------
2316       subroutine read_compar
2317 !
2318 ! Read molecular data
2319 !
2320       use conform_compar, only:alloc_compar_arrays
2321       use control_data, only:pdbref
2322       use geometry_data, only:deg2rad,rad2deg
2323 !      implicit none
2324 !      include 'DIMENSIONS'
2325 !      include 'DIMENSIONS.ZSCOPT'
2326 !      include 'DIMENSIONS.COMPAR'
2327 !      include 'DIMENSIONS.FREE'
2328 !      include 'COMMON.IOUNITS'
2329 !      include 'COMMON.TIME1'
2330 !      include 'COMMON.SBRIDGE'
2331 !      include 'COMMON.CONTROL'
2332 !      include 'COMMON.COMPAR'
2333 !      include 'COMMON.CHAIN'
2334 !      include 'COMMON.HEADER'
2335 !      include 'COMMON.GEO'
2336 !      include 'COMMON.FREE'
2337       character(len=320) :: controlcard !,ucase
2338       character(len=64) :: wfile
2339 !el      integer ilen
2340 !el      external ilen
2341       integer :: i,j,k
2342 !elwrite(iout,*)"jestesmy w read_compar"
2343       call card_concat(controlcard,.true.)
2344       pdbref=(index(controlcard,'PDBREF').gt.0)
2345       call reada(controlcard,'CUTOFF_UP',rmscut_base_up,4.0d0)
2346       call reada(controlcard,'CUTOFF_LOW',rmscut_base_low,3.0d0)
2347       call reada(controlcard,'RMSUP_LIM',rmsup_lim,4.0d0)
2348       call reada(controlcard,'RMSUPUP_LIM',rmsupup_lim,7.5d0)
2349       verbose = index(controlcard,"VERBOSE").gt.0
2350       lgrp=index(controlcard,"STATIN").gt.0
2351       lgrp_out=index(controlcard,"STATOUT").gt.0
2352       merge_helices=index(controlcard,"DONT_MERGE_HELICES").eq.0
2353       binary = index(controlcard,"BINARY").gt.0
2354       rmscut_base_up=rmscut_base_up/50
2355       rmscut_base_low=rmscut_base_low/50
2356       call reada(controlcard,"FRAC_SEC",frac_sec,0.66666666d0)
2357       call readi(controlcard,'NLEVEL',nlevel,1)
2358       if (nlevel.lt.0) then
2359         allocate(nfrag(2))
2360         call alloc_compar_arrays(maxfrag,1)
2361         goto 121
2362       else
2363         allocate(nfrag(nlevel))
2364       endif
2365 ! Read the data pertaining to elementary fragments (level 1)
2366       call readi(controlcard,'NFRAG',nfrag(1),0)
2367       write(iout,*)"nfrag(1)",nfrag(1)
2368       call alloc_compar_arrays(nfrag(1),nlevel)
2369       do j=1,nfrag(1)
2370         call card_concat(controlcard,.true.)
2371         write (iout,*) controlcard(:ilen(controlcard))
2372         call readi(controlcard,'NPIECE',npiece(j,1),0)
2373         call readi(controlcard,'N_SHIFT1',n_shift(1,j,1),0)
2374         call readi(controlcard,'N_SHIFT2',n_shift(2,j,1),0)
2375         call reada(controlcard,'ANGCUT',ang_cut(j),50.0d0)
2376         call reada(controlcard,'MAXANG',ang_cut1(j),360.0d0)
2377         call reada(controlcard,'FRAC_MIN',frac_min(j),0.666666d0)
2378         call reada(controlcard,'NC_FRAC',nc_fragm(j,1),0.5d0)
2379         call readi(controlcard,'NC_REQ',nc_req_setf(j,1),0)
2380         call readi(controlcard,'RMS',irms(j,1),0)
2381         call readi(controlcard,'LOCAL',iloc(j),1)
2382         call readi(controlcard,'ELCONT',ielecont(j,1),1)
2383         if (ielecont(j,1).eq.0) then
2384           call readi(controlcard,'SCCONT',isccont(j,1),1)
2385         endif
2386         ang_cut(j)=ang_cut(j)*deg2rad
2387         ang_cut1(j)=ang_cut1(j)*deg2rad
2388         do k=1,npiece(j,1)
2389           call card_concat(controlcard,.true.)
2390           call readi(controlcard,'IFRAG1',ifrag(1,k,j),0)
2391           call readi(controlcard,'IFRAG2',ifrag(2,k,j),0)
2392         enddo
2393         write(iout,*)"j",j," npiece",npiece(j,1)," ifrag",&
2394           (ifrag(1,k,j),ifrag(2,k,j),&
2395          k=1,npiece(j,1))," ang_cut",ang_cut(j)*rad2deg,&
2396           " ang_cut1",ang_cut1(j)*rad2deg
2397         write(iout,*)"n_shift",n_shift(1,j,1),n_shift(2,j,1)
2398         write(iout,*)"nc_frac",nc_fragm(j,1)," nc_req",nc_req_setf(j,1)
2399         write(iout,*)"irms",irms(j,1)," ielecont",ielecont(j,1),&
2400           " ilocal",iloc(j)," isccont",isccont(j,1)
2401       enddo
2402 ! Read data pertaning to higher levels
2403       do i=2,nlevel
2404         call card_concat(controlcard,.true.)
2405         call readi(controlcard,'NFRAG',NFRAG(i),0)
2406         write (iout,*) "i",i," nfrag",nfrag(i)
2407         do j=1,nfrag(i)
2408           call card_concat(controlcard,.true.)
2409           if (i.eq.2) then
2410             call readi(controlcard,'ELCONT',ielecont(j,i),0)
2411             if (ielecont(j,i).eq.0) then
2412               call readi(controlcard,'SCCONT',isccont(j,i),1)
2413             endif
2414             call readi(controlcard,'RMS',irms(j,i),0)
2415           else
2416             ielecont(j,i)=0
2417             isccont(j,i)=0
2418             irms(j,i)=1
2419           endif
2420           call readi(controlcard,'NPIECE',npiece(j,i),0)
2421           call readi(controlcard,'N_SHIFT1',n_shift(1,j,i),0)
2422           call readi(controlcard,'N_SHIFT2',n_shift(2,j,i),0)
2423           call multreadi(controlcard,'IPIECE',ipiece(1,j,i),&
2424             npiece(j,i),0)
2425           call reada(controlcard,'NC_FRAC',nc_fragm(j,i),0.5d0)
2426           call readi(controlcard,'NC_REQ',nc_req_setf(j,i),0)
2427           write(iout,*) "j",j," npiece",npiece(j,i)," n_shift",&
2428             n_shift(1,j,i),n_shift(2,j,i)," ielecont",ielecont(j,i),&
2429             " isccont",isccont(j,i)," irms",irms(j,i)
2430           write(iout,*) "ipiece",(ipiece(k,j,i),k=1,npiece(j,i))
2431           write(iout,*)"n_shift",n_shift(1,j,i),n_shift(2,j,i)
2432           write(iout,*)"nc_frac",nc_fragm(j,i),&
2433            " nc_req",nc_req_setf(j,i)
2434         enddo
2435       enddo
2436       if (binary) write (iout,*) "Classes written in binary format."
2437       return
2438   121 continue
2439       call reada(controlcard,'ANGCUT_HEL',angcut_hel,50.0d0)
2440       call reada(controlcard,'MAXANG_HEL',angcut1_hel,60.0d0)
2441       call reada(controlcard,'ANGCUT_BET',angcut_bet,90.0d0)
2442       call reada(controlcard,'MAXANG_BET',angcut1_bet,360.0d0)
2443       call reada(controlcard,'ANGCUT_STRAND',angcut_strand,90.0d0)
2444       call reada(controlcard,'MAXANG_STRAND',angcut1_strand,60.0d0)
2445       call reada(controlcard,'FRAC_MIN',frac_min_set,0.666666d0)
2446       call reada(controlcard,'NC_FRAC_HEL',ncfrac_hel,0.5d0)
2447       call readi(controlcard,'NC_REQ_HEL',ncreq_hel,0)
2448       call reada(controlcard,'NC_FRAC_BET',ncfrac_bet,0.5d0)
2449       call reada(controlcard,'NC_FRAC_PAIR',ncfrac_pair,0.3d0)
2450       call readi(controlcard,'NC_REQ_BET',ncreq_bet,0)
2451       call readi(controlcard,'NC_REQ_PAIR',ncreq_pair,0)
2452       call readi(controlcard,'NSHIFT_HEL',nshift_hel,3)
2453       call readi(controlcard,'NSHIFT_BET',nshift_bet,3)
2454       call readi(controlcard,'NSHIFT_STRAND',nshift_strand,3)
2455       call readi(controlcard,'NSHIFT_PAIR',nshift_pair,3)
2456       call readi(controlcard,'RMS_SINGLE',irms_single,0)
2457       call readi(controlcard,'CONT_SINGLE',icont_single,1)
2458       call readi(controlcard,'LOCAL_SINGLE',iloc_single,1)
2459       call readi(controlcard,'RMS_PAIR',irms_pair,0)
2460       call readi(controlcard,'CONT_PAIR',icont_pair,1)
2461       call readi(controlcard,'SPLIT_BET',isplit_bet,0)
2462       angcut_hel=angcut_hel*deg2rad
2463       angcut1_hel=angcut1_hel*deg2rad
2464       angcut_bet=angcut_bet*deg2rad
2465       angcut1_bet=angcut1_bet*deg2rad
2466       angcut_strand=angcut_strand*deg2rad
2467       angcut1_strand=angcut1_strand*deg2rad
2468       write (iout,*) "Automatic detection of structural elements"
2469       write (iout,*) 'NC_FRAC_HEL',ncfrac_hel,' NC_REQ_HEL',ncreq_hel,&
2470                      ' NC_FRAC_BET',ncfrac_bet,' NC_REQ_BET',ncreq_bet,&
2471                  ' RMS_SINGLE',irms_single,' CONT_SINGLE',icont_single,&
2472                  ' NC_FRAC_PAIR',ncfrac_pair,' NC_REQ_PAIR',ncreq_pair,&
2473         ' RMS_PAIR',irms_pair,' CONT_PAIR',icont_pair,&
2474         ' SPLIT_BET',isplit_bet
2475       write (iout,*) 'NSHIFT_HEL',nshift_hel,' NSHIFT_BET',nshift_bet,&
2476         ' NSHIFT_STRAND',nshift_strand,' NSHIFT_PAIR',nshift_pair
2477       write (iout,*) 'ANGCUT_HEL',angcut_hel*rad2deg,&
2478         ' MAXANG_HEL',angcut1_hel*rad2deg
2479       write (iout,*) 'ANGCUT_BET',angcut_bet*rad2deg,&
2480                      ' MAXANG_BET',angcut1_bet*rad2deg
2481       write (iout,*) 'ANGCUT_STRAND',angcut_strand*rad2deg,&
2482                      ' MAXANG_STRAND',angcut1_strand*rad2deg
2483       write (iout,*) 'FRAC_MIN',frac_min_set
2484       return
2485       end subroutine read_compar
2486 !--------------------------------------------------------------------------------
2487 ! read_ref_str.F
2488 !--------------------------------------------------------------------------------
2489       subroutine read_ref_structure(*)
2490 !
2491 ! Read the reference structure from the PDB file or from a PDB file or in the form of the dihedral
2492 ! angles.
2493 !
2494       use control_data, only:pdbref
2495       use geometry_data, only:nres,cref,c,dc,nsup,dc_norm,nend_sup,&
2496                               nstart_sup,nstart_seq,nperm,nres0
2497       use energy_data, only:nct,nnt,icont_ref,ncont_ref,itype
2498       use compare, only:seq_comp !,contact,elecont
2499       use geometry, only:chainbuild,dist
2500       use io_config, only:readpdb
2501 !
2502       use conform_compar, only:contact,elecont
2503 !      implicit none
2504 !      include 'DIMENSIONS'
2505 !      include 'DIMENSIONS.ZSCOPT'
2506 !      include 'DIMENSIONS.COMPAR'
2507 !      include 'COMMON.IOUNITS'
2508 !      include 'COMMON.GEO'
2509 !      include 'COMMON.VAR'
2510 !      include 'COMMON.INTERACT'
2511 !      include 'COMMON.LOCAL'
2512 !      include 'COMMON.NAMES'
2513 !      include 'COMMON.CHAIN'
2514 !      include 'COMMON.FFIELD'
2515 !      include 'COMMON.SBRIDGE'
2516 !      include 'COMMON.HEADER'
2517 !      include 'COMMON.CONTROL'
2518 !      include 'COMMON.CONTACTS1'
2519 !      include 'COMMON.PEPTCONT'
2520 !      include 'COMMON.TIME1'
2521 !      include 'COMMON.COMPAR'
2522       character(len=4) :: sequence(nres)
2523 !el      integer rescode
2524 !el      real(kind=8) :: x(maxvar)
2525       integer :: itype_pdb(nres)
2526 !el      logical seq_comp
2527       integer :: i,j,k,nres_pdb,iaux
2528       real(kind=8) :: ddsc !el,dist
2529       integer :: kkk !,ilen
2530 !el      external ilen
2531 !
2532       nres0=nres
2533       write (iout,*) "pdbref",pdbref
2534       if (pdbref) then
2535         read(inp,'(a)') pdbfile
2536         write (iout,'(2a,1h.)') 'PDB data will be read from file ',&
2537           pdbfile(:ilen(pdbfile))
2538         open(ipdbin,file=pdbfile,status='old',err=33)
2539         goto 34
2540   33    write (iout,'(a)') 'Error opening PDB file.'
2541         return 1
2542   34    continue
2543         do i=1,nres
2544           itype_pdb(i)=itype(i)
2545         enddo
2546
2547         call readpdb
2548
2549         do i=1,nres
2550           iaux=itype_pdb(i)
2551           itype_pdb(i)=itype(i)
2552           itype(i)=iaux
2553         enddo
2554         close (ipdbin)
2555         do kkk=1,nperm
2556         nres_pdb=nres
2557         nres=nres0
2558         nstart_seq=nnt
2559         if (nsup.le.(nct-nnt+1)) then
2560           do i=0,nct-nnt+1-nsup
2561             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),&
2562               nsup)) then
2563               do j=nnt+nsup-1,nnt,-1
2564                 do k=1,3
2565                   cref(k,nres+j+i,kkk)=cref(k,nres_pdb+j,kkk)
2566                 enddo
2567               enddo
2568               do j=nnt+nsup-1,nnt,-1
2569                 do k=1,3
2570                   cref(k,j+i,kkk)=cref(k,j,kkk)
2571                 enddo
2572                 phi_ref(j+i)=phi_ref(j)
2573                 theta_ref(j+i)=theta_ref(j)
2574                 alph_ref(j+i)=alph_ref(j)
2575                 omeg_ref(j+i)=omeg_ref(j)
2576               enddo
2577 #ifdef DEBUG
2578               do j=nnt,nct
2579                 write (iout,'(i5,3f10.5,5x,3f10.5)') &
2580                   j,(cref(k,j,kkk),k=1,3),(cref(k,j+nres,kkk),k=1,3)
2581               enddo
2582 #endif
2583               nstart_seq=nnt+i
2584               nstart_sup=nnt+i
2585               goto 111
2586             endif
2587           enddo
2588           write (iout,'(a)') &
2589                   'Error - sequences to be superposed do not match.'
2590           return 1
2591         else
2592           do i=0,nsup-(nct-nnt+1)
2593             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),&
2594               nct-nnt+1)) &
2595             then
2596               nstart_sup=nstart_sup+i
2597               nsup=nct-nnt+1
2598               goto 111
2599             endif
2600           enddo
2601           write (iout,'(a)') &
2602                   'Error - sequences to be superposed do not match.'
2603         endif
2604         enddo
2605   111   continue
2606         write (iout,'(a,i5)') &
2607          'Experimental structure begins at residue',nstart_seq
2608       else
2609         call read_angles(inp,*38)
2610         goto 39
2611    38   write (iout,'(a)') 'Error reading reference structure.'
2612         return 1
2613    39   call chainbuild
2614         kkk=1
2615         nstart_sup=nnt
2616         nstart_seq=nnt
2617         nsup=nct-nnt+1
2618         do i=1,2*nres
2619           do j=1,3
2620             cref(j,i,kkk)=c(j,i)
2621           enddo
2622         enddo
2623       endif
2624       nend_sup=nstart_sup+nsup-1
2625       do i=1,2*nres
2626         do j=1,3
2627           c(j,i)=cref(j,i,kkk)
2628         enddo
2629       enddo
2630       do i=1,nres
2631         do j=1,3
2632           dc(j,nres+i)=cref(j,nres+i,kkk)-cref(j,i,kkk)
2633         enddo
2634         if (itype(i).ne.10) then
2635           ddsc = dist(i,nres+i)
2636           do j=1,3
2637             dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
2638           enddo
2639         else
2640           do j=1,3
2641             dc_norm(j,nres+i)=0.0d0
2642           enddo
2643         endif
2644 !        write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
2645 !         " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
2646 !         dc_norm(3,nres+i)**2
2647         do j=1,3
2648           dc(j,i)=c(j,i+1)-c(j,i)
2649         enddo
2650         ddsc = dist(i,i+1)
2651         do j=1,3
2652           dc_norm(j,i)=dc(j,i)/ddsc
2653         enddo
2654       enddo
2655 !      print *,"Calling contact"
2656       call contact(.true.,ncont_ref,icont_ref(1,1),&
2657         nstart_sup,nend_sup)
2658 !      print *,"Calling elecont"
2659       call elecont(.true.,ncont_pept_ref,&
2660          icont_pept_ref(1,1),&
2661          nstart_sup,nend_sup)
2662        write (iout,'(a,i3,a,i3,a,i3,a)') &
2663           'Number of residues to be superposed:',nsup,&
2664           ' (from residue',nstart_sup,' to residue',&
2665           nend_sup,').'
2666       return
2667       end subroutine read_ref_structure
2668 !--------------------------------------------------------------------------------
2669 ! geomout.F
2670 !--------------------------------------------------------------------------------
2671       subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev)
2672
2673       use geometry_data, only:nres,c
2674       use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype
2675 !      implicit real*8 (a-h,o-z)
2676 !      include 'DIMENSIONS'
2677 !      include 'DIMENSIONS.ZSCOPT'
2678 !      include 'COMMON.CHAIN'
2679 !      include 'COMMON.INTERACT'
2680 !      include 'COMMON.NAMES'
2681 !      include 'COMMON.IOUNITS'
2682 !      include 'COMMON.HEADER'
2683 !      include 'COMMON.SBRIDGE'
2684       character(len=50) :: tytul
2685       character(len=1),dimension(10) :: chainid=reshape((/'A','B','C',&
2686                       'D','E','F','G','H','I','J'/),shape(chainid))
2687       integer,dimension(nres) :: ica !(maxres)
2688       real(kind=8) :: temp,efree,etot,entropy,rmsdev
2689       integer :: ii,i,j,iti,ires,iatom,ichain
2690       write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')&
2691         ii,temp,rmsdev
2692       write (ipdb,'("REMARK DIMENSIONLESS FREE ENERGY",1pe15.5)') &
2693         efree
2694       write (ipdb,'("REMARK ENERGY",1pe15.5," ENTROPY",1pe15.5)') &
2695         etot,entropy
2696       iatom=0
2697       ichain=1
2698       ires=0
2699       do i=nnt,nct
2700         iti=itype(i)
2701         if (iti.eq.ntyp1) then
2702           ichain=ichain+1
2703           ires=0
2704           write (ipdb,'(a)') 'TER'
2705         else
2706         ires=ires+1
2707         iatom=iatom+1
2708         ica(i)=iatom
2709         write (ipdb,10) iatom,restyp(iti),chainid(ichain),&
2710            ires,(c(j,i),j=1,3)
2711         if (iti.ne.10) then
2712           iatom=iatom+1
2713           write (ipdb,20) iatom,restyp(iti),chainid(ichain),&
2714             ires,(c(j,nres+i),j=1,3)
2715         endif
2716         endif
2717       enddo
2718       write (ipdb,'(a)') 'TER'
2719       do i=nnt,nct-1
2720         if (itype(i).eq.ntyp1) cycle
2721         if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
2722           write (ipdb,30) ica(i),ica(i+1)
2723         else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
2724           write (ipdb,30) ica(i),ica(i+1),ica(i)+1
2725         else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
2726           write (ipdb,30) ica(i),ica(i)+1
2727         endif
2728       enddo
2729       if (itype(nct).ne.10) then
2730         write (ipdb,30) ica(nct),ica(nct)+1
2731       endif
2732       do i=1,nss
2733         write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
2734       enddo
2735       write (ipdb,'(a6)') 'ENDMDL'
2736   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
2737   20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
2738   30  FORMAT ('CONECT',8I5)
2739       return
2740       end subroutine pdboutW
2741 #endif
2742 !------------------------------------------------------------------------------
2743       end module io_wham
2744 !-----------------------------------------------------------------------------
2745 !-----------------------------------------------------------------------------
2746