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