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