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('BONDPAR_NUCL',bondname_nucl)
  55       open (ibond_nucl,file=bondname_nucl,status='old')
  56       call mygetenv('THETPAR',thetname)
  57       open (ithep,file=thetname,status='old')
  58       call mygetenv('ROTPAR',rotname)
  59       open (irotam,file=rotname,status='old')
  60       call mygetenv('TORPAR',torname)
  61       open (itorp,file=torname,status='old')
  62       call mygetenv('TORDPAR',tordname)
  63       open (itordp,file=tordname,status='old')
  64       call mygetenv('FOURIER',fouriername)
  65       open (ifourier,file=fouriername,status='old')
  66       call mygetenv('SCCORPAR',sccorname)
  67       open (isccor,file=sccorname,status='old')
  68       call mygetenv('ELEPAR',elename)
  69       open (ielep,file=elename,status='old')
  70       call mygetenv('SIDEPAR',sidename)
  71       open (isidep,file=sidename,status='old')
  72       call mygetenv('SIDEP',sidepname)
  73       open (isidep1,file=sidepname,status="old")
  74       call mygetenv('THETPAR_NUCL',thetname_nucl)
  75       open (ithep_nucl,file=thetname_nucl,status='old')
  76       call mygetenv('ROTPAR_NUCL',rotname_nucl)
  77       open (irotam_nucl,file=rotname_nucl,status='old')
  78       call mygetenv('TORPAR_NUCL',torname_nucl)
  79       open (itorp_nucl,file=torname_nucl,status='old')
  80       call mygetenv('TORDPAR_NUCL',tordname_nucl)
  81       open (itordp_nucl,file=tordname_nucl,status='old')
  82       call mygetenv('SIDEPAR_NUCL',sidename_nucl)
  83       open (isidep_nucl,file=sidename_nucl,status='old')
  84       call mygetenv('SIDEPAR_SCBASE',sidename_scbase)
  85       open (isidep_scbase,file=sidename_scbase,status='old')
  86       call mygetenv('PEPPAR_PEPBASE',pepname_pepbase)
  87       open (isidep_pepbase,file=pepname_pepbase,status='old')
  88       call mygetenv('SCPAR_PHOSPH',pepname_scpho)
  89       open (isidep_scpho,file=pepname_scpho,status='old')
  90       call mygetenv('PEPPAR_PHOSPH',pepname_peppho)
  91       open (isidep_peppho,file=pepname_peppho,status='old')
  92       call mygetenv('LIPTRANPAR',liptranname)
  93       open (iliptranpar,file=liptranname,status='old',action='read')
  94       call mygetenv('TUBEPAR',tubename)
  95       open (itube,file=tubename,status='old',action='read')
  96       call mygetenv('IONPAR',ionname)
  97       open (iion,file=ionname,status='old',action='read')
  98
  99       call mygetenv('SCPPAR_NUCL',scpname_nucl)
 100       open (iscpp_nucl,file=scpname_nucl,status='old')
 101
 102
 103 #ifndef OLDSCP
 104 !
 105 ! 8/9/01 In the newest version SCp interaction constants are read from a file
 106 ! Use -DOLDSCP to use hard-coded constants instead.
 107 !
 108       call mygetenv('SCPPAR',scpname)
 109       open (iscpp,file=scpname,status='old')
 110 #endif
 111 #ifdef MPL
 112       if (MyID.eq.BossID) then
 113       MyRank = MyID/fgProcs
 114 #endif
 115 #ifdef MPI
 116       print *,'OpenUnits: processor',MyRank
 117       call numstr(MyRank,liczba)
 118       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//liczba
 119 #else
 120       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
 121 #endif
 122       open(iout,file=outname,status='unknown')
 123       write (iout,'(80(1h-))')
 124       write (iout,'(30x,a)') "FILE ASSIGNMENT"
 125       write (iout,'(80(1h-))')
 126       write (iout,*) "Input file                      : ",&
 127         prefix(:ilen(prefix))//'.inp'
 128       write (iout,*) "Output file                     : ",&
 129         outname(:ilen(outname))
 130       write (iout,*)
 131       write (iout,*) "Sidechain potential file        : ",&
 132         sidename(:ilen(sidename))
 133 #ifndef OLDSCP
 134       write (iout,*) "SCp potential file              : ",&
 135         scpname(:ilen(scpname))
 136 #endif
 137       write (iout,*) "Electrostatic potential file    : ",&
 138         elename(:ilen(elename))
 139       write (iout,*) "Cumulant coefficient file       : ",&
 140         fouriername(:ilen(fouriername))
 141       write (iout,*) "Torsional parameter file        : ",&
 142         torname(:ilen(torname))
 143       write (iout,*) "Double torsional parameter file : ",&
 144         tordname(:ilen(tordname))
 145       write (iout,*) "Backbone-rotamer parameter file : ",&
 146         sccorname(:ilen(sccorname))
 147       write (iout,*) "Bond & inertia constant file    : ",&
 148         bondname(:ilen(bondname))
 149       write (iout,*) "Bending parameter file          : ",&
 150         thetname(:ilen(thetname))
 151       write (iout,*) "Rotamer parameter file          : ",&
 152         rotname(:ilen(rotname))
 153       write (iout,'(80(1h-))')
 154       write (iout,*)
 155       return
 156       end subroutine openunits
 157 !-----------------------------------------------------------------------------
 158 ! molread_zs.F
 159 !-----------------------------------------------------------------------------
 160       subroutine molread(*)
 161 !
 162 ! Read molecular data.
 163 !
 164       use energy_data
 165       use geometry_data, only:nres,deg2rad,c,dc,nres_molec
 166       use control_data, only:iscode
 167       use io_base, only:rescode
 168       use control, only:setup_var,init_int_table
 169       use geometry, only:alloc_geo_arrays
 170       use energy, only:alloc_ener_arrays
 171 !      implicit real*8 (a-h,o-z)
 172 !      include 'DIMENSIONS'
 173 !      include 'DIMENSIONS.ZSCOPT'
 174 !      include 'COMMON.IOUNITS'
 175 !      include 'COMMON.GEO'
 176 !      include 'COMMON.VAR'
 177 !      include 'COMMON.INTERACT'
 178 !      include 'COMMON.LOCAL'
 179 !      include 'COMMON.NAMES'
 180 !      include 'COMMON.CHAIN'
 181 !      include 'COMMON.FFIELD'
 182 !      include 'COMMON.SBRIDGE'
 183 !      include 'COMMON.TORCNSTR'
 184 !      include 'COMMON.CONTROL'
 185       character(len=4),dimension(:),allocatable :: sequence !(nres)
 186 !el      integer :: rescode
 187 !el      real(kind=8) :: x(maxvar)
 188       character(len=320) :: controlcard !,ucase
 189       integer,dimension(nres,5) :: itype_pdb !(maxres)
 190       integer :: i,j,i1,i2,it1,it2,mnum
 191       real(kind=8) :: scalscp
 192 !el      logical :: seq_comp
 193       write(iout,*) "START MOLREAD"
 194       call flush(iout)
 195       do i=1,5
 196        nres_molec(i)=0
 197        print *,"nres_molec, initial",nres_molec(i)
 198       enddo
 199
 200       call card_concat(controlcard,.true.)
 201       call reada(controlcard,'SCAL14',scal14,0.4d0)
 202       call reada(controlcard,'SCALSCP',scalscp,1.0d0)
 203       call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
 204       call reada(controlcard,'TEMP0',temp0,300.0d0) !el
 205       call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
 206       r0_corr=cutoff_corr-delt_corr
 207       call readi(controlcard,"NRES",nres_molec(1),0)
 208       call readi(controlcard,"NRES_NUCL",nres_molec(2),0)
 209       call readi(controlcard,"NRES_CAT",nres_molec(5),0)
 210       nres=0
 211       do i=1,5
 212        nres=nres_molec(i)+nres
 213       enddo
 214 !      allocate(sequence(nres+1))
 215 !el znamy juz ilosc reszt wiec mozna zaalokowac tablice do liczenia enerii
 216       call alloc_geo_arrays
 217       call alloc_ener_arrays
 218 ! alokacja dodatkowych tablic, ktore w unresie byly alokowanie w locie
 219 !----------------------------
 220       allocate(c(3,2*nres+2))
 221       allocate(dc(3,0:2*nres+2))
 222       allocate(itype(nres+2,5))
 223       allocate(itel(nres+2))
 224       if (.not. allocated(molnum)) allocate(molnum(nres+2))
 225 !
 226 ! Zero out tableis.
 227       do i=1,2*nres+2
 228         do j=1,3
 229           c(j,i)=0.0D0
 230           dc(j,i)=0.0D0
 231         enddo
 232       enddo
 233       do i=1,nres+2
 234         molnum(i)=0
 235         do j=1,5
 236         itype(i,j)=0
 237         enddo
 238         itel(i)=0
 239       enddo
 240 !--------------------------
 241 !
 242       allocate(sequence(nres+1))
 243       iscode=index(controlcard,"ONE_LETTER")
 244       if (nres.le.0) then
 245         write (iout,*) "Error: no residues in molecule"
 246         return 1
 247       endif
 248       if (nres.gt.maxres) then
 249         write (iout,*) "Error: too many residues",nres,maxres
 250       endif
 251       write(iout,*) 'nres=',nres
 252 ! HERE F**** CHANGE
 253 ! Read sequence of the protein
 254       if (iscode.gt.0) then
 255         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 256       else
 257         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 258       endif
 259 ! Convert sequence to numeric code
 260       do i=1,nres_molec(1)
 261         mnum=1
 262         molnum(i)=1
 263         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
 264       enddo
 265       do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
 266         mnum=2
 267         molnum(i)=2
 268         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
 269       enddo
 270       do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
 271         mnum=5
 272         molnum(i)=5
 273         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
 274       enddo
 275
 276       write (iout,*) "Numeric code:"
 277       write (iout,'(20i4)') (itype(i,molnum(i)),i=1,nres)
 278       do i=1,nres-1
 279         mnum=molnum(i)
 280 #ifdef PROCOR
 281         if (itype(i,mnum).eq.ntyp1_molec(mnum) .or. itype(i+1,mnum).eq.ntyp1_molec(mnum)) then
 282 #else
 283         if (itype(i,mnum).eq.ntyp1_molec(mnum)) then
 284 #endif
 285           itel(i)=0
 286 #ifdef PROCOR
 287         else if (iabs(itype(i+1,mnum)).ne.20) then
 288 #else
 289         else if (iabs(itype(i,mnum)).ne.20) then
 290 #endif
 291           itel(i)=1
 292         else
 293           itel(i)=2
 294         endif
 295       enddo
 296        write (iout,*) "ITEL"
 297        do i=1,nres-1
 298          mnum=molnum(i)
 299          write (iout,*) i,itype(i,mnum),itel(i)
 300        enddo
 301       write(iout,*)
 302       call read_bridge
 303
 304       if (with_dihed_constr) then
 305
 306       read (inp,*) ndih_constr
 307       if (ndih_constr.gt.0) then
 308         read (inp,*) ftors
 309         write (iout,*) 'FTORS',ftors
 310         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 311         write (iout,*) &
 312          'There are',ndih_constr,' constraints on phi angles.'
 313         do i=1,ndih_constr
 314           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 315         enddo
 316         do i=1,ndih_constr
 317           phi0(i)=deg2rad*phi0(i)
 318           drange(i)=deg2rad*drange(i)
 319         enddo
 320       endif
 321
 322       endif
 323
 324       nnt=1
 325       nct=nres
 326       if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
 327       if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
 328       write(iout,*) 'NNT=',NNT,' NCT=',NCT
 329       call setup_var
 330       call init_int_table
 331       if (ns.gt.0) then
 332         write (iout,'(/a,i3,a)') 'The chain contains',ns,&
 333         ' disulfide-bridging cysteines.'
 334         write (iout,'(20i4)') (iss(i),i=1,ns)
 335         write (iout,'(/a/)') 'Pre-formed links are:'
 336         do i=1,nss
 337           mnum=molnum(i)
 338           i1=ihpb(i)-nres
 339           i2=jhpb(i)-nres
 340           it1=itype(i1,mnum)
 341           it2=itype(i2,mnum)
 342          write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
 343           restyp(it1,molnum(i1)),'(',i1,') -- ',restyp(it2,molnum(i2)),'(',i2,')',&
 344           dhpb(i),ebr,forcon(i)
 345         enddo
 346       endif
 347       write (iout,'(a)')
 348       return
 349       end subroutine molread
 350 !-----------------------------------------------------------------------------
 351 ! parmread.F
 352 !-----------------------------------------------------------------------------
 353       subroutine parmread(iparm,*)
 354 #else
 355       subroutine parmread
 356 #endif
 357 !
 358 ! Read the parameters of the probability distributions of the virtual-bond
 359 ! valence angles and the side chains and energy parameters.
 360 !
 361       use wham_data
 362
 363       use geometry_data
 364       use energy_data
 365       use control_data, only: maxterm,maxlor,maxterm_sccor,& !maxtor
 366           maxtermd_1,maxtermd_2,tor_mode,scelemode !,maxthetyp,maxthetyp1
 367       use MD_data
 368 !el      use MPI_data
 369 !el      use map_data
 370       use io_config, only: printmat
 371       use control, only: getenv_loc
 372
 373 #ifdef MPI
 374       use MPI_data
 375       include "mpif.h"
 376       integer :: IERROR
 377 #endif
 378 !      implicit real*8 (a-h,o-z)
 379 !      include 'DIMENSIONS'
 380 !      include 'DIMENSIONS.ZSCOPT'
 381 !      include 'DIMENSIONS.FREE'
 382 !      include 'COMMON.IOUNITS'
 383 !      include 'COMMON.CHAIN'
 384 !      include 'COMMON.INTERACT'
 385 !      include 'COMMON.GEO'
 386 !      include 'COMMON.LOCAL'
 387 !      include 'COMMON.TORSION'
 388 !      include 'COMMON.FFIELD'
 389 !      include 'COMMON.NAMES'
 390 !      include 'COMMON.SBRIDGE'
 391 !      include 'COMMON.WEIGHTS'
 392 !      include 'COMMON.ENEPS'
 393 !      include 'COMMON.SCCOR'
 394 !      include 'COMMON.SCROT'
 395 !      include 'COMMON.FREE'
 396       character(len=1) :: t1,t2,t3
 397       character(len=1) :: onelett(4) = (/"G","A","P","D"/)
 398       character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
 399       logical :: lprint,SPLIT_FOURIERTOR
 400       real(kind=8),dimension(3,3,maxlob) :: blower      !(3,3,maxlob)
 401       character(len=800) :: controlcard
 402       character(len=256) :: bondname_t,thetname_t,rotname_t,torname_t,&
 403         tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,&
 404         sccorname_t
 405 !el      integer ilen
 406 !el   external ilen
 407       character(len=16) :: key
 408       integer :: iparm,nkcctyp
 409 !el      real(kind=8) :: ip,mp
 410       real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,epsijlip,&
 411                 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,    &
 412                 epspeptube,epssctube,sigmapeptube,      &
 413                 sigmasctube
 414
 415       real(kind=8) :: v0ij,v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,rjunk,&
 416                 res1
 417       integer :: i,j,ichir1,ichir2,k,l,m,kk,ii,mm,junk,lll,ll,llll,n,jj
 418       integer :: nlobi,iblock,maxinter,iscprol,ncatprotparm
 419       character*3 string
 420
 421 !
 422 ! Body
 423 !
 424 ! Set LPRINT=.TRUE. for debugging
 425       dwa16=2.0d0**(1.0d0/6.0d0)
 426       lprint=.false.
 427       itypro=20
 428 ! Assign virtual-bond length
 429       vbl=3.8D0
 430       vblinv=1.0D0/vbl
 431       vblinv2=vblinv*vblinv
 432 #ifndef CLUSTER
 433       call card_concat(controlcard,.true.)
 434       wname(4)="WCORRH"
 435 !el
 436 allocate(ww(max_eneW))
 437       do i=1,n_eneW
 438         key = wname(i)(:ilen(wname(i)))
 439         call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
 440       enddo
 441
 442       write (iout,*) "iparm",iparm," myparm",myparm
 443 ! If reading not own parameters, skip assignment
 444
 445       if (iparm.eq.myparm .or. .not.separate_parset) then
 446
 447 !
 448 ! Setup weights for UNRES
 449 !
 450       wsc=ww(1)
 451       wscp=ww(2)
 452       welec=ww(3)
 453       wcorr=ww(4)
 454       wcorr5=ww(5)
 455       wcorr6=ww(6)
 456       wel_loc=ww(7)
 457       wturn3=ww(8)
 458       wturn4=ww(9)
 459       wturn6=ww(10)
 460       wang=ww(11)
 461       wscloc=ww(12)
 462       wtor=ww(13)
 463       wtor_d=ww(14)
 464       wvdwpp=ww(16)
 465       wbond=ww(18)
 466       wsccor=ww(19)
 467       wcatcat=ww(42)
 468       wcatprot=ww(41)
 469
 470       endif
 471 !
 472 !el------
 473       allocate(weights(n_ene))
 474       weights(1)=wsc
 475       weights(2)=wscp
 476       weights(3)=welec
 477       weights(4)=wcorr
 478       weights(5)=wcorr5
 479       weights(6)=wcorr6
 480       weights(7)=wel_loc
 481       weights(8)=wturn3
 482       weights(9)=wturn4
 483       weights(10)=wturn6
 484       weights(11)=wang
 485       weights(12)=wscloc
 486       weights(13)=wtor
 487       weights(14)=wtor_d
 488       weights(15)=0 !wstrain !
 489       weights(16)=0 !wvdwpp !
 490       weights(17)=wbond
 491       weights(18)=0 !scal14 !
 492       weights(21)=wsccor
 493       weights(42)=wcatprot
 494       weights(41)=wcatcat
 495
 496 ! el--------
 497       call card_concat(controlcard,.false.)
 498
 499 ! Return if not own parameters
 500
 501       if (iparm.ne.myparm .and. separate_parset) return
 502
 503       call reads(controlcard,"BONDPAR",bondname_t,bondname)
 504       open (ibond,file=bondname_t,status='old')
 505       rewind(ibond)
 506       call reads(controlcard,"THETPAR",thetname_t,thetname)
 507       open (ithep,file=thetname_t,status='old')
 508       rewind(ithep)
 509       call reads(controlcard,"ROTPAR",rotname_t,rotname)
 510       open (irotam,file=rotname_t,status='old')
 511       rewind(irotam)
 512       call reads(controlcard,"TORPAR",torname_t,torname)
 513       open (itorp,file=torname_t,status='old')
 514       rewind(itorp)
 515       call reads(controlcard,"TORDPAR",tordname_t,tordname)
 516       open (itordp,file=tordname_t,status='old')
 517       rewind(itordp)
 518       call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
 519       open (isccor,file=sccorname_t,status='old')
 520       rewind(isccor)
 521       call reads(controlcard,"FOURIER",fouriername_t,fouriername)
 522       open (ifourier,file=fouriername_t,status='old')
 523       rewind(ifourier)
 524       call reads(controlcard,"ELEPAR",elename_t,elename)
 525       open (ielep,file=elename_t,status='old')
 526       rewind(ielep)
 527       call reads(controlcard,"SIDEPAR",sidename_t,sidename)
 528       open (isidep,file=sidename_t,status='old')
 529       rewind(isidep)
 530       call reads(controlcard,"SCPPAR",scpname_t,scpname)
 531       open (iscpp,file=scpname_t,status='old')
 532       rewind(iscpp)
 533       call getenv_loc('IONPAR',ionname)
 534       open (iion,file=ionname,status='old')
 535       rewind(iion)
 536       write (iout,*) "Parameter set:",iparm
 537       write (iout,*) "Energy-term weights:"
 538       do i=1,n_eneW
 539         write (iout,'(a16,f10.5)') wname(i),ww(i)
 540       enddo
 541       write (iout,*) "Sidechain potential file        : ",&
 542         sidename_t(:ilen(sidename_t))
 543 #ifndef OLDSCP
 544       write (iout,*) "SCp potential file              : ",&
 545         scpname_t(:ilen(scpname_t))
 546 #endif
 547       write (iout,*) "Electrostatic potential file    : ",&
 548         elename_t(:ilen(elename_t))
 549       write (iout,*) "Cumulant coefficient file       : ",&
 550         fouriername_t(:ilen(fouriername_t))
 551       write (iout,*) "Torsional parameter file        : ",&
 552         torname_t(:ilen(torname_t))
 553       write (iout,*) "Double torsional parameter file : ",&
 554         tordname_t(:ilen(tordname_t))
 555       write (iout,*) "Backbone-rotamer parameter file : ",&
 556         sccorname(:ilen(sccorname))
 557       write (iout,*) "Bond & inertia constant file    : ",&
 558         bondname_t(:ilen(bondname_t))
 559       write (iout,*) "Bending parameter file          : ",&
 560         thetname_t(:ilen(thetname_t))
 561       write (iout,*) "Rotamer parameter file          : ",&
 562         rotname_t(:ilen(rotname_t))
 563 #endif
 564 !
 565 ! Read the virtual-bond parameters, masses, and moments of inertia
 566 ! and Stokes' radii of the peptide group and side chains
 567 !
 568       allocate(dsc(ntyp1)) !(ntyp1)
 569       allocate(dsc_inv(ntyp1)) !(ntyp1)
 570       allocate(nbondterm(ntyp)) !(ntyp)
 571       allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 572       allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 573       allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
 574       allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
 575       allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
 576
 577 !el      allocate(msc(ntyp+1)) !(ntyp+1)
 578 !el      allocate(isc(ntyp+1)) !(ntyp+1)
 579 !el      allocate(restok(ntyp+1)) !(ntyp+1)
 580       allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 581
 582 #ifdef CRYST_BOND
 583       read (ibond,*) vbldp0,akp
 584       do i=1,ntyp
 585         nbondterm(i)=1
 586         read (ibond,*) vbldsc0(1,i),aksc(1,i)
 587         dsc(i) = vbldsc0(1,i)
 588         if (i.eq.10) then
 589           dsc_inv(i)=0.0D0
 590         else
 591           dsc_inv(i)=1.0D0/dsc(i)
 592         endif
 593       enddo
 594 #else
 595       read (ibond,*) ijunk,vbldp0,vbldpDUM,akp,rjunk
 596       do i=1,ntyp
 597         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
 598          j=1,nbondterm(i))
 599         dsc(i) = vbldsc0(1,i)
 600         if (i.eq.10) then
 601           dsc_inv(i)=0.0D0
 602         else
 603           dsc_inv(i)=1.0D0/dsc(i)
 604         endif
 605       enddo
 606 #endif
 607       if (lprint) then
 608         write(iout,'(/a/)')"Force constants virtual bonds:"
 609         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',&
 610          'inertia','Pstok'
 611         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
 612         do i=1,ntyp
 613           write (iout,'(a10,i3,6f10.5)') restyp(i,molnum(i)),nbondterm(i),&
 614             vbldsc0(1,i),aksc(1,i),abond0(1,i)
 615           do j=2,nbondterm(i)
 616             write (iout,'(13x,3f10.5)') &
 617               vbldsc0(j,i),aksc(j,i),abond0(j,i)
 618           enddo
 619         enddo
 620       endif
 621             if (.not. allocated(msc)) allocate(msc(ntyp1,5))
 622             if (.not. allocated(restok)) allocate(restok(ntyp1,5))
 623
 624             do i=1,ntyp_molec(5)
 625              read(iion,*) msc(i,5),restok(i,5)
 626              print *,msc(i,5),restok(i,5)
 627             enddo
 628             ip(5)=0.2
 629 !            isc(5)=0.2
 630             read (iion,*) ncatprotparm
 631             allocate(catprm(ncatprotparm,4))
 632             do k=1,4
 633             read (iion,*)  (catprm(i,k),i=1,ncatprotparm)
 634             enddo
 635             print *, catprm
 636       read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
 637       do i=1,ntyp_molec(2)
 638         nbondterm_nucl(i)=1
 639         read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),rjunk,restok(i,2)
 640 !        dsc(i) = vbldsc0_nucl(1,i)
 641 !        if (i.eq.10) then
 642 !          dsc_inv(i)=0.0D0
 643 !        else
 644 !          dsc_inv(i)=1.0D0/dsc(i)
 645 !        endif
 646       enddo
 647
 648
 649 !----------------------------------------------------
 650       allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
 651       allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp))      !(-ntyp:ntyp)
 652       allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
 653       allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
 654       allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
 655       allocate(gthet(3,-ntyp:ntyp))     !(3,-ntyp:ntyp)
 656       do i=-ntyp,ntyp
 657         a0thet(i)=0.0D0
 658         do j=1,2
 659          do ichir1=-1,1
 660           do ichir2=-1,1
 661           athet(j,i,ichir1,ichir2)=0.0D0
 662           bthet(j,i,ichir1,ichir2)=0.0D0
 663           enddo
 664          enddo
 665         enddo
 666         do j=0,3
 667           polthet(j,i)=0.0D0
 668         enddo
 669         do j=1,3
 670           gthet(j,i)=0.0D0
 671         enddo
 672         theta0(i)=0.0D0
 673         sig0(i)=0.0D0
 674         sigc0(i)=0.0D0
 675       enddo
 676 !elwrite(iout,*) "parmread kontrol"
 677
 678 #ifdef CRYST_THETA
 679 !
 680 ! Read the parameters of the probability distribution/energy expression
 681 ! of the virtual-bond valence angles theta
 682 !
 683       do i=1,ntyp
 684         read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),&
 685           (bthet(j,i,1,1),j=1,2)
 686         read (ithep,*) (polthet(j,i),j=0,3)
 687 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 688         read (ithep,*) (gthet(j,i),j=1,3)
 689 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 690         read (ithep,*) theta0(i),sig0(i),sigc0(i)
 691         sigc0(i)=sigc0(i)**2
 692 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 693       enddo
 694 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 695       do i=1,ntyp
 696       athet(1,i,1,-1)=athet(1,i,1,1)
 697       athet(2,i,1,-1)=athet(2,i,1,1)
 698       bthet(1,i,1,-1)=-bthet(1,i,1,1)
 699       bthet(2,i,1,-1)=-bthet(2,i,1,1)
 700       athet(1,i,-1,1)=-athet(1,i,1,1)
 701       athet(2,i,-1,1)=-athet(2,i,1,1)
 702       bthet(1,i,-1,1)=bthet(1,i,1,1)
 703       bthet(2,i,-1,1)=bthet(2,i,1,1)
 704       enddo
 705 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 706       do i=-ntyp,-1
 707       a0thet(i)=a0thet(-i)
 708       athet(1,i,-1,-1)=athet(1,-i,1,1)
 709       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 710       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 711       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 712       athet(1,i,-1,1)=athet(1,-i,1,1)
 713       athet(2,i,-1,1)=-athet(2,-i,1,1)
 714       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 715       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 716       athet(1,i,1,-1)=-athet(1,-i,1,1)
 717       athet(2,i,1,-1)=athet(2,-i,1,1)
 718       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 719       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 720       theta0(i)=theta0(-i)
 721       sig0(i)=sig0(-i)
 722       sigc0(i)=sigc0(-i)
 723        do j=0,3
 724         polthet(j,i)=polthet(j,-i)
 725        enddo
 726        do j=1,3
 727          gthet(j,i)=gthet(j,-i)
 728        enddo
 729       enddo
 730 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 731       close (ithep)
 732 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 733       if (lprint) then
 734 !       write (iout,'(a)')
 735 !    &    'Parameters of the virtual-bond valence angles:'
 736 !       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
 737 !    & '    ATHETA0   ','         A1   ','        A2    ',
 738 !    & '        B1    ','         B2   '
 739 !       do i=1,ntyp
 740 !         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 741 !    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
 742 !       enddo
 743 !       write (iout,'(/a/9x,5a/79(1h-))')
 744 !    & 'Parameters of the expression for sigma(theta_c):',
 745 !    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
 746 !    & '     ALPH3    ','    SIGMA0C   '
 747 !       do i=1,ntyp
 748 !         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 749 !    &      (polthet(j,i),j=0,3),sigc0(i)
 750 !       enddo
 751 !       write (iout,'(/a/9x,5a/79(1h-))')
 752 !    & 'Parameters of the second gaussian:',
 753 !    & '    THETA0    ','     SIGMA0   ','        G1    ',
 754 !    & '        G2    ','         G3   '
 755 !       do i=1,ntyp
 756 !         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
 757 !    &       sig0(i),(gthet(j,i),j=1,3)
 758 !       enddo
 759         write (iout,'(a)') &
 760           'Parameters of the virtual-bond valence angles:'
 761         write (iout,'(/a/9x,5a/79(1h-))') &
 762        'Coefficients of expansion',&
 763        '     theta0   ','    a1*10^2   ','   a2*10^2    ',&
 764        '   b1*10^1    ','    b2*10^1   '
 765         do i=1,ntyp
 766           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
 767               a0thet(i),(100*athet(j,i,1,1),j=1,2),&
 768               (10*bthet(j,i,1,1),j=1,2)
 769         enddo
 770         write (iout,'(/a/9x,5a/79(1h-))') &
 771        'Parameters of the expression for sigma(theta_c):',&
 772        ' alpha0       ','  alph1       ',' alph2        ',&
 773        ' alhp3        ','   sigma0c    '
 774         do i=1,ntyp
 775           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
 776             (polthet(j,i),j=0,3),sigc0(i)
 777         enddo
 778         write (iout,'(/a/9x,5a/79(1h-))') &
 779        'Parameters of the second gaussian:',&
 780        '    theta0    ','  sigma0*10^2 ','      G1*10^-1',&
 781        '        G2    ','   G3*10^1    '
 782         do i=1,ntyp
 783           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
 784              100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 785         enddo
 786       endif
 787 #else
 788 !
 789 ! Read the parameters of Utheta determined from ab initio surfaces
 790 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 791 !
 792       allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
 793 !      write (iout,*) "tu dochodze"
 794       IF (tor_mode.eq.0) THEN
 795 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 796       read (ithep,*) nthetyp,ntheterm,ntheterm2,&
 797         ntheterm3,nsingle,ndouble
 798       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 799
 800 !----------------------------------------------------
 801 !      allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
 802       allocate(aa0thet(-nthetyp-1:nthetyp+1,&
 803         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 804 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 805       allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
 806         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 807 !(maxtheterm,-maxthetyp1:maxthetyp1,&
 808 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 809       allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 810         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 811       allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 812         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 813       allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 814         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 815       allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 816         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 817 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
 818 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 819       allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
 820         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 821       allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
 822         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 823 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
 824 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 825
 826
 827       read (ithep,*) (ithetyp(i),i=1,ntyp1)
 828       do i=-ntyp1,-1
 829         ithetyp(i)=-ithetyp(-i)
 830       enddo
 831 !      write (iout,*) "tu dochodze"
 832       aa0thet(:,:,:,:)=0.0d0
 833       aathet(:,:,:,:,:)=0.0d0
 834       bbthet(:,:,:,:,:,:)=0.0d0
 835       ccthet(:,:,:,:,:,:)=0.0d0
 836       ddthet(:,:,:,:,:,:)=0.0d0
 837       eethet(:,:,:,:,:,:)=0.0d0
 838       ffthet(:,:,:,:,:,:,:)=0.0d0
 839       ggthet(:,:,:,:,:,:,:)=0.0d0
 840
 841       do iblock=1,2
 842       do i=0,nthetyp
 843         do j=-nthetyp,nthetyp
 844           do k=-nthetyp,nthetyp
 845             read (ithep,'(6a)') res1
 846             read (ithep,*) aa0thet(i,j,k,iblock)
 847             read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
 848             read (ithep,*) &
 849              ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
 850               (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
 851               (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
 852               (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
 853               ll=1,ntheterm2)
 854             read (ithep,*) &
 855             (((ffthet(llll,lll,ll,i,j,k,iblock),&
 856                ffthet(lll,llll,ll,i,j,k,iblock),&
 857                ggthet(llll,lll,ll,i,j,k,iblock),&
 858                ggthet(lll,llll,ll,i,j,k,iblock),&
 859                llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 860           enddo
 861         enddo
 862       enddo
 863 !
 864 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
 865 ! coefficients of theta-and-gamma-dependent terms are zero.
 866 !
 867       do i=1,nthetyp
 868         do j=1,nthetyp
 869           do l=1,ntheterm
 870             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
 871             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
 872           enddo
 873           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
 874           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
 875         enddo
 876         do l=1,ntheterm
 877           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 878         enddo
 879         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
 880       enddo
 881       enddo
 882 ! Substitution for D aminoacids from symmetry.
 883       do iblock=1,2
 884       do i=-nthetyp,0
 885         do j=-nthetyp,nthetyp
 886           do k=-nthetyp,nthetyp
 887            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
 888            do l=1,ntheterm
 889            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
 890            enddo
 891            do ll=1,ntheterm2
 892             do lll=1,nsingle
 893             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
 894             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
 895             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
 896             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
 897             enddo
 898           enddo
 899           do ll=1,ntheterm3
 900            do lll=2,ndouble
 901             do llll=1,lll-1
 902             ffthet(llll,lll,ll,i,j,k,iblock)= &
 903             ffthet(llll,lll,ll,-i,-j,-k,iblock)
 904             ffthet(lll,llll,ll,i,j,k,iblock)= &
 905             ffthet(lll,llll,ll,-i,-j,-k,iblock)
 906             ggthet(llll,lll,ll,i,j,k,iblock)= &
 907             -ggthet(llll,lll,ll,-i,-j,-k,iblock)
 908             ggthet(lll,llll,ll,i,j,k,iblock)= &
 909             -ggthet(lll,llll,ll,-i,-j,-k,iblock)
 910             enddo !ll
 911            enddo  !lll
 912           enddo   !llll
 913          enddo    !k
 914         enddo     !j
 915        enddo      !i
 916       enddo       !iblock
 917
 918 !
 919 ! Control printout of the coefficients of virtual-bond-angle potentials
 920 !
 921       do iblock=1,2
 922       if (lprint) then
 923         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
 924         do i=1,nthetyp+1
 925           do j=1,nthetyp+1
 926             do k=1,nthetyp+1
 927               write (iout,'(//4a)') &
 928                'Type ',onelett(i),onelett(j),onelett(k)
 929               write (iout,'(//a,10x,a)') " l","a[l]"
 930               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
 931               write (iout,'(i2,1pe15.5)') &
 932                  (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
 933             do l=1,ntheterm2
 934               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))') &
 935                 "b",l,"c",l,"d",l,"e",l
 936               do m=1,nsingle
 937                 write (iout,'(i2,4(1pe15.5))') m,&
 938                 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
 939                 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
 940               enddo
 941             enddo
 942             do l=1,ntheterm3
 943               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))') &
 944                 "f+",l,"f-",l,"g+",l,"g-",l
 945               do m=2,ndouble
 946                 do n=1,m-1
 947                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
 948                     ffthet(n,m,l,i,j,k,iblock),&
 949                     ffthet(m,n,l,i,j,k,iblock),&
 950                     ggthet(n,m,l,i,j,k,iblock),&
 951                     ggthet(m,n,l,i,j,k,iblock)
 952                 enddo
 953               enddo
 954             enddo
 955           enddo
 956         enddo
 957       enddo
 958       call flush(iout)
 959       endif
 960       enddo
 961       ELSE
 962 !C here will be the apropriate recalibrating for D-aminoacid
 963       read (ithep,*) nthetyp
 964       allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
 965       allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
 966       do i=-nthetyp+1,nthetyp-1
 967         read (ithep,*) nbend_kcc_Tb(i)
 968         do j=0,nbend_kcc_Tb(i)
 969           read (ithep,*) ijunk,v1bend_chyb(j,i)
 970         enddo
 971       enddo
 972       if (lprint) then
 973         write (iout,'(a)') &
 974          "Parameters of the valence-only potentials"
 975         do i=-nthetyp+1,nthetyp-1
 976         write (iout,'(2a)') "Type ",toronelet(i)
 977         do k=0,nbend_kcc_Tb(i)
 978           write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
 979         enddo
 980         enddo
 981       endif
 982       ENDIF ! TOR_MODE
 983
 984 #endif
 985 !--------------- Reading theta parameters for nucleic acid-------
 986       read (ithep_nucl,*) nthetyp_nucl,ntheterm_nucl,&
 987       ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
 988       nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
 989       allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
 990       allocate(aa0thet_nucl(nthetyp_nucl+1,&
 991         nthetyp_nucl+1,nthetyp_nucl+1))
 992 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 993       allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
 994         nthetyp_nucl+1,nthetyp_nucl+1))
 995 !(maxtheterm,-maxthetyp1:maxthetyp1,&
 996 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 997       allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
 998         nthetyp_nucl+1,nthetyp_nucl+1))
 999       allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1000         nthetyp_nucl+1,nthetyp_nucl+1))
1001       allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1002         nthetyp_nucl+1,nthetyp_nucl+1))
1003       allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1004         nthetyp_nucl+1,nthetyp_nucl+1))
1005 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1006 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1007       allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1008          nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1009       allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1010          nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1011
1012 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1013 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1014
1015       read (ithep_nucl,*) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1016
1017       aa0thet_nucl(:,:,:)=0.0d0
1018       aathet_nucl(:,:,:,:)=0.0d0
1019       bbthet_nucl(:,:,:,:,:)=0.0d0
1020       ccthet_nucl(:,:,:,:,:)=0.0d0
1021       ddthet_nucl(:,:,:,:,:)=0.0d0
1022       eethet_nucl(:,:,:,:,:)=0.0d0
1023       ffthet_nucl(:,:,:,:,:,:)=0.0d0
1024       ggthet_nucl(:,:,:,:,:,:)=0.0d0
1025
1026       do i=1,nthetyp_nucl
1027         do j=1,nthetyp_nucl
1028           do k=1,nthetyp_nucl
1029             read (ithep_nucl,'(3a)') t1,t2,t3
1030             read (ithep_nucl,*) aa0thet_nucl(i,j,k)
1031             read (ithep_nucl,*)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1032             read (ithep_nucl,*) &
1033             (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1034             (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1035             (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1036             (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1037             read (ithep_nucl,*) &
1038            (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1039               ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1040               llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1041           enddo
1042         enddo
1043       enddo
1044
1045
1046 !-------------------------------------------
1047       allocate(nlob(ntyp1)) !(ntyp1)
1048       allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1049       allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1050       allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1051
1052       do i=1,ntyp
1053         do j=1,maxlob
1054           bsc(j,i)=0.0D0
1055           nlob(i)=0
1056         enddo
1057       enddo
1058       nlob(ntyp1)=0
1059       dsc(ntyp1)=0.0D0
1060
1061       do i=-ntyp,ntyp
1062         do j=1,maxlob
1063           do k=1,3
1064             censc(k,j,i)=0.0D0
1065           enddo
1066           do k=1,3
1067             do l=1,3
1068               gaussc(l,k,j,i)=0.0D0
1069             enddo
1070           enddo
1071         enddo
1072       enddo
1073
1074 #ifdef CRYST_SC
1075 !
1076 ! Read the parameters of the probability distribution/energy expression
1077 ! of the side chains.
1078 !
1079       do i=1,ntyp
1080 !c      write (iout,*) "tu dochodze",i
1081         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
1082         if (i.eq.10) then
1083           dsc_inv(i)=0.0D0
1084         else
1085           dsc_inv(i)=1.0D0/dsc(i)
1086         endif
1087         if (i.ne.10) then
1088         do j=1,nlob(i)
1089           do k=1,3
1090             do l=1,3
1091               blower(l,k,j)=0.0D0
1092             enddo
1093           enddo
1094         enddo
1095         bsc(1,i)=0.0D0
1096         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
1097         censc(1,1,-i)=censc(1,1,i)
1098         censc(2,1,-i)=censc(2,1,i)
1099         censc(3,1,-i)=-censc(3,1,i)
1100         do j=2,nlob(i)
1101           read (irotam,*) bsc(j,i)
1102           read (irotam,*) (censc(k,j,i),k=1,3),&
1103                                        ((blower(k,l,j),l=1,k),k=1,3)
1104         censc(1,j,-i)=censc(1,j,i)
1105         censc(2,j,-i)=censc(2,j,i)
1106         censc(3,j,-i)=-censc(3,j,i)
1107 ! BSC is amplitude of Gaussian
1108         enddo
1109         do j=1,nlob(i)
1110           do k=1,3
1111             do l=1,k
1112               akl=0.0D0
1113               do m=1,3
1114                 akl=akl+blower(k,m,j)*blower(l,m,j)
1115               enddo
1116               gaussc(k,l,j,i)=akl
1117               gaussc(l,k,j,i)=akl
1118              if (((k.eq.3).and.(l.ne.3)) &
1119               .or.((l.eq.3).and.(k.ne.3))) then
1120                 gaussc(k,l,j,-i)=-akl
1121                 gaussc(l,k,j,-i)=-akl
1122               else
1123                 gaussc(k,l,j,-i)=akl
1124                 gaussc(l,k,j,-i)=akl
1125               endif
1126             enddo
1127           enddo
1128         enddo
1129         endif
1130       enddo
1131       close (irotam)
1132       if (lprint) then
1133         write (iout,'(/a)') 'Parameters of side-chain local geometry'
1134         do i=1,ntyp
1135           nlobi=nlob(i)
1136           if (nlobi.gt.0) then
1137           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
1138            ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1139 !          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1140 !          write (iout,'(a,f10.4,4(16x,f10.4))')
1141 !     &                             'Center  ',(bsc(j,i),j=1,nlobi)
1142 !          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
1143            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1144                                    'log h',(bsc(j,i),j=1,nlobi)
1145            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1146           'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1147 !          write (iout,'(a)')
1148 !         do j=1,nlobi
1149 !           ind=0
1150 !           do k=1,3
1151 !             do l=1,k
1152 !              ind=ind+1
1153 !              blower(k,l,j)=gaussc(ind,j,i)
1154 !             enddo
1155 !           enddo
1156 !         enddo
1157           do k=1,3
1158             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1159                        ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1160           enddo
1161           endif
1162         enddo
1163       endif
1164 #else
1165 !
1166 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1167 ! added by Urszula Kozlowska 07/11/2007
1168 !
1169       allocate(sc_parmin(65,ntyp))      !(maxsccoef,ntyp)
1170
1171       do i=1,ntyp
1172         read (irotam,*)
1173        if (i.eq.10) then
1174          read (irotam,*)
1175        else
1176          do j=1,65
1177            read(irotam,*) sc_parmin(j,i)
1178          enddo
1179        endif
1180       enddo
1181 #endif
1182       close(irotam)
1183 !---------reading nucleic acid parameters for rotamers-------------------
1184       allocate(sc_parmin_nucl(9,ntyp_molec(2)))      !(maxsccoef,ntyp)
1185       do i=1,ntyp_molec(2)
1186         read (irotam_nucl,*)
1187         do j=1,9
1188           read(irotam_nucl,*) sc_parmin_nucl(j,i)
1189         enddo
1190       enddo
1191       close(irotam_nucl)
1192       if (lprint) then
1193         write (iout,*)
1194         write (iout,*) "Base rotamer parameters"
1195         do i=1,ntyp_molec(2)
1196           write (iout,'(a)') restyp(i,2)
1197           write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1198         enddo
1199       endif
1200
1201
1202       read (ifourier,*) nloctyp
1203 !el write(iout,*)"nloctyp",nloctyp
1204       SPLIT_FOURIERTOR = nloctyp.lt.0
1205       nloctyp = iabs(nloctyp)
1206 #ifdef NEWCORR
1207       if (.not.allocated(itype2loc)) allocate(itype2loc(-ntyp1:ntyp1))
1208        print *,"shape",shape(itype2loc)
1209       allocate(iloctyp(-nloctyp:nloctyp))
1210       allocate(bnew1(3,2,-nloctyp:nloctyp))
1211       allocate(bnew2(3,2,-nloctyp:nloctyp))
1212       allocate(ccnew(3,2,-nloctyp:nloctyp))
1213       allocate(ddnew(3,2,-nloctyp:nloctyp))
1214       allocate(e0new(3,-nloctyp:nloctyp))
1215       allocate(eenew(2,2,2,-nloctyp:nloctyp))
1216       allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1217       allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1218       allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1219       allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1220       allocate(e0newtor(3,-nloctyp:nloctyp))
1221       allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1222
1223       read (ifourier,*) (itype2loc(i),i=1,ntyp)
1224       read (ifourier,*) (iloctyp(i),i=0,nloctyp-1)
1225       itype2loc(ntyp1)=nloctyp
1226       iloctyp(nloctyp)=ntyp1
1227       do i=1,ntyp1
1228         itype2loc(-i)=-itype2loc(i)
1229       enddo
1230 #else
1231       allocate(iloctyp(-nloctyp:nloctyp))
1232       allocate(itype2loc(-ntyp1:ntyp1))
1233       iloctyp(0)=10
1234       iloctyp(1)=9
1235       iloctyp(2)=20
1236       iloctyp(3)=ntyp1
1237 #endif
1238       do i=1,nloctyp
1239         iloctyp(-i)=-iloctyp(i)
1240       enddo
1241 !c      write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1242 !c      write (iout,*) "nloctyp",nloctyp,
1243 !c     &  " iloctyp",(iloctyp(i),i=0,nloctyp)
1244 !c      write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1245 !c      write (iout,*) "nloctyp",nloctyp,
1246 !c     &  " iloctyp",(iloctyp(i),i=0,nloctyp)
1247 #ifdef NEWCORR
1248       do i=0,nloctyp-1
1249 !c             write (iout,*) "NEWCORR",i
1250         read (ifourier,*)
1251         do ii=1,3
1252           do j=1,2
1253             read (ifourier,*) bnew1(ii,j,i)
1254           enddo
1255         enddo
1256 !c             write (iout,*) "NEWCORR BNEW1"
1257 !c             write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1258         do ii=1,3
1259           do j=1,2
1260             read (ifourier,*) bnew2(ii,j,i)
1261           enddo
1262         enddo
1263 !c             write (iout,*) "NEWCORR BNEW2"
1264 !c             write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1265         do kk=1,3
1266           read (ifourier,*) ccnew(kk,1,i)
1267           read (ifourier,*) ccnew(kk,2,i)
1268         enddo
1269 !c             write (iout,*) "NEWCORR CCNEW"
1270 !c             write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1271         do kk=1,3
1272           read (ifourier,*) ddnew(kk,1,i)
1273           read (ifourier,*) ddnew(kk,2,i)
1274         enddo
1275 !c             write (iout,*) "NEWCORR DDNEW"
1276 !c             write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1277         do ii=1,2
1278           do jj=1,2
1279             do kk=1,2
1280               read (ifourier,*) eenew(ii,jj,kk,i)
1281             enddo
1282           enddo
1283         enddo
1284 !c             write (iout,*) "NEWCORR EENEW1"
1285 !c             write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1286         do ii=1,3
1287           read (ifourier,*) e0new(ii,i)
1288         enddo
1289 !c             write (iout,*) (e0new(ii,i),ii=1,3)
1290       enddo
1291 !c             write (iout,*) "NEWCORR EENEW"
1292       do i=0,nloctyp-1
1293         do ii=1,3
1294           ccnew(ii,1,i)=ccnew(ii,1,i)/2
1295           ccnew(ii,2,i)=ccnew(ii,2,i)/2
1296           ddnew(ii,1,i)=ddnew(ii,1,i)/2
1297           ddnew(ii,2,i)=ddnew(ii,2,i)/2
1298         enddo
1299       enddo
1300       do i=1,nloctyp-1
1301         do ii=1,3
1302           bnew1(ii,1,-i)= bnew1(ii,1,i)
1303           bnew1(ii,2,-i)=-bnew1(ii,2,i)
1304           bnew2(ii,1,-i)= bnew2(ii,1,i)
1305           bnew2(ii,2,-i)=-bnew2(ii,2,i)
1306         enddo
1307         do ii=1,3
1308 !c          ccnew(ii,1,i)=ccnew(ii,1,i)/2
1309 !c          ccnew(ii,2,i)=ccnew(ii,2,i)/2
1310 !c          ddnew(ii,1,i)=ddnew(ii,1,i)/2
1311 !c          ddnew(ii,2,i)=ddnew(ii,2,i)/2
1312           ccnew(ii,1,-i)=ccnew(ii,1,i)
1313           ccnew(ii,2,-i)=-ccnew(ii,2,i)
1314           ddnew(ii,1,-i)=ddnew(ii,1,i)
1315           ddnew(ii,2,-i)=-ddnew(ii,2,i)
1316         enddo
1317         e0new(1,-i)= e0new(1,i)
1318         e0new(2,-i)=-e0new(2,i)
1319         e0new(3,-i)=-e0new(3,i)
1320         do kk=1,2
1321           eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1322           eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1323           eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1324           eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1325         enddo
1326       enddo
1327       if (lprint) then
1328         write (iout,'(a)') "Coefficients of the multibody terms"
1329         do i=-nloctyp+1,nloctyp-1
1330           write (iout,*) "Type: ",onelet(iloctyp(i))
1331           write (iout,*) "Coefficients of the expansion of B1"
1332           do j=1,2
1333             write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1334           enddo
1335           write (iout,*) "Coefficients of the expansion of B2"
1336           do j=1,2
1337             write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1338           enddo
1339           write (iout,*) "Coefficients of the expansion of C"
1340           write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1341           write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1342           write (iout,*) "Coefficients of the expansion of D"
1343           write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1344           write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1345           write (iout,*) "Coefficients of the expansion of E"
1346           write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1347           do j=1,2
1348             do k=1,2
1349               write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1350             enddo
1351           enddo
1352         enddo
1353       endif
1354       IF (SPLIT_FOURIERTOR) THEN
1355       do i=0,nloctyp-1
1356 !c             write (iout,*) "NEWCORR TOR",i
1357         read (ifourier,*)
1358         do ii=1,3
1359           do j=1,2
1360             read (ifourier,*) bnew1tor(ii,j,i)
1361           enddo
1362         enddo
1363 !c             write (iout,*) "NEWCORR BNEW1 TOR"
1364 !c             write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1365         do ii=1,3
1366           do j=1,2
1367             read (ifourier,*) bnew2tor(ii,j,i)
1368           enddo
1369         enddo
1370 !c             write (iout,*) "NEWCORR BNEW2 TOR"
1371 !c             write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1372         do kk=1,3
1373           read (ifourier,*) ccnewtor(kk,1,i)
1374           read (ifourier,*) ccnewtor(kk,2,i)
1375         enddo
1376 !c             write (iout,*) "NEWCORR CCNEW TOR"
1377 !c             write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1378         do kk=1,3
1379           read (ifourier,*) ddnewtor(kk,1,i)
1380           read (ifourier,*) ddnewtor(kk,2,i)
1381         enddo
1382 !c             write (iout,*) "NEWCORR DDNEW TOR"
1383 !c             write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1384         do ii=1,2
1385           do jj=1,2
1386             do kk=1,2
1387               read (ifourier,*) eenewtor(ii,jj,kk,i)
1388             enddo
1389           enddo
1390         enddo
1391 !c         write (iout,*) "NEWCORR EENEW1 TOR"
1392 !c         write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1393         do ii=1,3
1394           read (ifourier,*) e0newtor(ii,i)
1395         enddo
1396 !c             write (iout,*) (e0newtor(ii,i),ii=1,3)
1397       enddo
1398 !c             write (iout,*) "NEWCORR EENEW TOR"
1399       do i=0,nloctyp-1
1400         do ii=1,3
1401           ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1402           ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1403           ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1404           ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1405         enddo
1406       enddo
1407       do i=1,nloctyp-1
1408         do ii=1,3
1409           bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1410           bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1411           bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1412           bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1413         enddo
1414         do ii=1,3
1415           ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1416           ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1417           ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1418           ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1419         enddo
1420         e0newtor(1,-i)= e0newtor(1,i)
1421         e0newtor(2,-i)=-e0newtor(2,i)
1422         e0newtor(3,-i)=-e0newtor(3,i)
1423         do kk=1,2
1424           eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1425           eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1426           eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1427           eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1428         enddo
1429       enddo
1430       if (lprint) then
1431         write (iout,'(a)') &
1432          "Single-body coefficients of the torsional potentials"
1433         do i=-nloctyp+1,nloctyp-1
1434           write (iout,*) "Type: ",onelet(iloctyp(i))
1435           write (iout,*) "Coefficients of the expansion of B1tor"
1436           do j=1,2
1437             write (iout,'(3hB1(,i1,1h),3f10.5)') &
1438              j,(bnew1tor(k,j,i),k=1,3)
1439           enddo
1440           write (iout,*) "Coefficients of the expansion of B2tor"
1441           do j=1,2
1442             write (iout,'(3hB2(,i1,1h),3f10.5)') &
1443              j,(bnew2tor(k,j,i),k=1,3)
1444           enddo
1445           write (iout,*) "Coefficients of the expansion of Ctor"
1446           write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1447           write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1448           write (iout,*) "Coefficients of the expansion of Dtor"
1449           write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1450           write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1451           write (iout,*) "Coefficients of the expansion of Etor"
1452           write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1453           do j=1,2
1454             do k=1,2
1455               write (iout,'(1hE,2i1,2f10.5)') &
1456                j,k,(eenewtor(l,j,k,i),l=1,2)
1457             enddo
1458           enddo
1459         enddo
1460       endif
1461       ELSE
1462       do i=-nloctyp+1,nloctyp-1
1463         do ii=1,3
1464           do j=1,2
1465             bnew1tor(ii,j,i)=bnew1(ii,j,i)
1466           enddo
1467         enddo
1468         do ii=1,3
1469           do j=1,2
1470             bnew2tor(ii,j,i)=bnew2(ii,j,i)
1471           enddo
1472         enddo
1473         do ii=1,3
1474           ccnewtor(ii,1,i)=ccnew(ii,1,i)
1475           ccnewtor(ii,2,i)=ccnew(ii,2,i)
1476           ddnewtor(ii,1,i)=ddnew(ii,1,i)
1477           ddnewtor(ii,2,i)=ddnew(ii,2,i)
1478         enddo
1479       enddo
1480       ENDIF !SPLIT_FOURIER_TOR
1481 #else
1482       allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1483       allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1484       allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1485       allocate(b(13,-nloctyp-1:nloctyp+1))
1486       if (lprint) &
1487        write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1488       do i=0,nloctyp-1
1489         read (ifourier,*)
1490         read (ifourier,*) (b(ii,i),ii=1,13)
1491         if (lprint) then
1492         write (iout,*) 'Type ',onelet(iloctyp(i))
1493         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1494         endif
1495         if (i.gt.0) then
1496         b(2,-i)= b(2,i)
1497         b(3,-i)= b(3,i)
1498         b(4,-i)=-b(4,i)
1499         b(5,-i)=-b(5,i)
1500         endif
1501         CCold(1,1,i)= b(7,i)
1502         CCold(2,2,i)=-b(7,i)
1503         CCold(2,1,i)= b(9,i)
1504         CCold(1,2,i)= b(9,i)
1505         CCold(1,1,-i)= b(7,i)
1506         CCold(2,2,-i)=-b(7,i)
1507         CCold(2,1,-i)=-b(9,i)
1508         CCold(1,2,-i)=-b(9,i)
1509         DDold(1,1,i)= b(6,i)
1510         DDold(2,2,i)=-b(6,i)
1511         DDold(2,1,i)= b(8,i)
1512         DDold(1,2,i)= b(8,i)
1513         DDold(1,1,-i)= b(6,i)
1514         DDold(2,2,-i)=-b(6,i)
1515         DDold(2,1,-i)=-b(8,i)
1516         DDold(1,2,-i)=-b(8,i)
1517         EEold(1,1,i)= b(10,i)+b(11,i)
1518         EEold(2,2,i)=-b(10,i)+b(11,i)
1519         EEold(2,1,i)= b(12,i)-b(13,i)
1520         EEold(1,2,i)= b(12,i)+b(13,i)
1521         EEold(1,1,-i)= b(10,i)+b(11,i)
1522         EEold(2,2,-i)=-b(10,i)+b(11,i)
1523         EEold(2,1,-i)=-b(12,i)+b(13,i)
1524         EEold(1,2,-i)=-b(12,i)-b(13,i)
1525         write(iout,*) "TU DOCHODZE"
1526         print *,"JESTEM"
1527       enddo
1528       if (lprint) then
1529       write (iout,*)
1530       write (iout,*) &
1531       "Coefficients of the cumulants (independent of valence angles)"
1532       do i=-nloctyp+1,nloctyp-1
1533         write (iout,*) 'Type ',onelet(iloctyp(i))
1534         write (iout,*) 'B1'
1535         write(iout,'(2f10.5)') B(3,i),B(5,i)
1536         write (iout,*) 'B2'
1537         write(iout,'(2f10.5)') B(2,i),B(4,i)
1538         write (iout,*) 'CC'
1539         do j=1,2
1540           write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1541         enddo
1542         write(iout,*) 'DD'
1543         do j=1,2
1544           write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1545         enddo
1546         write(iout,*) 'EE'
1547         do j=1,2
1548           write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1549         enddo
1550       enddo
1551       endif
1552 #endif
1553 #ifdef CRYST_TOR
1554 !
1555 ! Read torsional parameters in old format
1556 !
1557       allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1558
1559       read (itorp,*) ntortyp,nterm_old
1560       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1561       read (itorp,*) (itortyp(i),i=1,ntyp)
1562
1563 !el from energy module--------
1564       allocate(v1(nterm_old,ntortyp,ntortyp))
1565       allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1566 !el---------------------------
1567
1568       do i=1,ntortyp
1569         do j=1,ntortyp
1570           read (itorp,'(a)')
1571           do k=1,nterm_old
1572             read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
1573           enddo
1574         enddo
1575       enddo
1576       close (itorp)
1577       if (lprint) then
1578         write (iout,'(/a/)') 'Torsional constants:'
1579         do i=1,ntortyp
1580           do j=1,ntortyp
1581             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1582             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1583           enddo
1584         enddo
1585       endif
1586
1587
1588 #else
1589 !
1590 ! Read torsional parameters
1591 !
1592       IF (TOR_MODE.eq.0) THEN
1593
1594       allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1595
1596       read (itorp,*) ntortyp
1597       read (itorp,*) (itortyp(i),i=1,ntyp)
1598       write (iout,*) 'ntortyp',ntortyp
1599
1600 !el from energy module---------
1601       allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1602       allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1603
1604       allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1605       allocate(vlor2(maxlor,ntortyp,ntortyp))
1606       allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1607       allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1608
1609       allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1610       allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1611 !el---------------------------
1612       do iblock=1,2
1613         do i=-ntortyp,ntortyp
1614           do j=-ntortyp,ntortyp
1615             nterm(i,j,iblock)=0
1616             nlor(i,j,iblock)=0
1617           enddo
1618         enddo
1619       enddo
1620 !el---------------------------
1621
1622       do iblock=1,2
1623       do i=-ntyp,-1
1624        itortyp(i)=-itortyp(-i)
1625       enddo
1626 !      write (iout,*) 'ntortyp',ntortyp
1627       do i=0,ntortyp-1
1628         do j=-ntortyp+1,ntortyp-1
1629           read (itorp,*) nterm(i,j,iblock),&
1630                 nlor(i,j,iblock)
1631           nterm(-i,-j,iblock)=nterm(i,j,iblock)
1632           nlor(-i,-j,iblock)=nlor(i,j,iblock)
1633           v0ij=0.0d0
1634           si=-1.0d0
1635           do k=1,nterm(i,j,iblock)
1636             read (itorp,*) kk,v1(k,i,j,iblock),&
1637             v2(k,i,j,iblock)
1638             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1639             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1640             v0ij=v0ij+si*v1(k,i,j,iblock)
1641             si=-si
1642          enddo
1643           do k=1,nlor(i,j,iblock)
1644             read (itorp,*) kk,vlor1(k,i,j),&
1645               vlor2(k,i,j),vlor3(k,i,j)
1646             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1647           enddo
1648           v0(i,j,iblock)=v0ij
1649           v0(-i,-j,iblock)=v0ij
1650         enddo
1651       enddo
1652       enddo
1653       close (itorp)
1654       if (lprint) then
1655         do iblock=1,2 !el
1656         write (iout,'(/a/)') 'Torsional constants:'
1657         do i=1,ntortyp
1658           do j=1,ntortyp
1659             write (iout,*) 'ityp',i,' jtyp',j
1660             write (iout,*) 'Fourier constants'
1661             do k=1,nterm(i,j,iblock)
1662               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1663               v2(k,i,j,iblock)
1664             enddo
1665             write (iout,*) 'Lorenz constants'
1666             do k=1,nlor(i,j,iblock)
1667               write (iout,'(3(1pe15.5))') &
1668                vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1669             enddo
1670           enddo
1671         enddo
1672         enddo
1673       endif
1674 !
1675 ! 6/23/01 Read parameters for double torsionals
1676 !
1677 !el from energy module------------
1678       allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1679       allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1680 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1681       allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1682       allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1683         !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1684       allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1685       allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1686         !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1687 !---------------------------------
1688
1689       do iblock=1,2
1690       do i=0,ntortyp-1
1691         do j=-ntortyp+1,ntortyp-1
1692           do k=-ntortyp+1,ntortyp-1
1693             read (itordp,'(3a1)') t1,t2,t3
1694 !              write (iout,*) "OK onelett",
1695 !     &         i,j,k,t1,t2,t3
1696
1697             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1698               .or. t3.ne.toronelet(k)) then
1699               write (iout,*) "Error in double torsional parameter file",&
1700                i,j,k,t1,t2,t3
1701 #ifdef MPI
1702               call MPI_Finalize(Ierror)
1703 #endif
1704                stop "Error in double torsional parameter file"
1705             endif
1706           read (itordp,*) ntermd_1(i,j,k,iblock),&
1707                ntermd_2(i,j,k,iblock)
1708             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1709             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1710             read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,&
1711                ntermd_1(i,j,k,iblock))
1712             read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,&
1713                ntermd_1(i,j,k,iblock))
1714             read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,&
1715                ntermd_1(i,j,k,iblock))
1716             read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,&
1717                ntermd_1(i,j,k,iblock))
1718 ! Martix of D parameters for one dimesional foureir series
1719             do l=1,ntermd_1(i,j,k,iblock)
1720              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1721              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1722              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1723              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1724 !            write(iout,*) "whcodze" ,
1725 !     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1726             enddo
1727             read (itordp,*) ((v2c(l,m,i,j,k,iblock),&
1728                v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1729                v2s(m,l,i,j,k,iblock),&
1730                m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1731 ! Martix of D parameters for two dimesional fourier series
1732             do l=1,ntermd_2(i,j,k,iblock)
1733              do m=1,l-1
1734              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1735              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1736              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1737              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1738              enddo!m
1739             enddo!l
1740           enddo!k
1741         enddo!j
1742       enddo!i
1743       enddo!iblock
1744       if (lprint) then
1745       write (iout,*)
1746       write (iout,*) 'Constants for double torsionals'
1747       do iblock=1,2
1748       do i=0,ntortyp-1
1749         do j=-ntortyp+1,ntortyp-1
1750           do k=-ntortyp+1,ntortyp-1
1751             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1752               ' nsingle',ntermd_1(i,j,k,iblock),&
1753               ' ndouble',ntermd_2(i,j,k,iblock)
1754             write (iout,*)
1755             write (iout,*) 'Single angles:'
1756             do l=1,ntermd_1(i,j,k,iblock)
1757               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1758                  v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1759                  v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1760                  v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1761             enddo
1762             write (iout,*)
1763             write (iout,*) 'Pairs of angles:'
1764             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1765             do l=1,ntermd_2(i,j,k,iblock)
1766               write (iout,'(i5,20f10.5)') &
1767                l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1768             enddo
1769             write (iout,*)
1770            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1771             do l=1,ntermd_2(i,j,k,iblock)
1772               write (iout,'(i5,20f10.5)') &
1773                l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1774                (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1775             enddo
1776             write (iout,*)
1777           enddo
1778         enddo
1779       enddo
1780       enddo
1781       endif
1782 #ifndef NEWCORR
1783       do i=1,ntyp1
1784         itype2loc(i)=itortyp(i)
1785       enddo
1786 #endif
1787       ELSE IF (TOR_MODE.eq.1) THEN
1788
1789 !C read valence-torsional parameters
1790       read (itorp,*) ntortyp
1791       nkcctyp=ntortyp
1792       write (iout,*) "Valence-torsional parameters read in ntortyp",&
1793         ntortyp
1794       read (itorp,*) (itortyp(i),i=1,ntyp)
1795       write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1796 #ifndef NEWCORR
1797       do i=1,ntyp1
1798         itype2loc(i)=itortyp(i)
1799       enddo
1800 #endif
1801       do i=-ntyp,-1
1802         itortyp(i)=-itortyp(-i)
1803       enddo
1804       do i=-ntortyp+1,ntortyp-1
1805         do j=-ntortyp+1,ntortyp-1
1806 !C first we read the cos and sin gamma parameters
1807           read (itorp,'(13x,a)') string
1808           write (iout,*) i,j,string
1809           read (itorp,*) &
1810          nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1811 !C           read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1812           do k=1,nterm_kcc(j,i)
1813             do l=1,nterm_kcc_Tb(j,i)
1814               do ll=1,nterm_kcc_Tb(j,i)
1815               read (itorp,*) ii,jj,kk, &
1816                v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1817               enddo
1818             enddo
1819           enddo
1820         enddo
1821       enddo
1822       ELSE
1823 #ifdef NEWCORR
1824 !c AL 4/8/16: Calculate coefficients from one-body parameters
1825       ntortyp=nloctyp
1826       allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1827       allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
1828       allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
1829       allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
1830       allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
1831
1832       do i=-ntyp1,ntyp1
1833        print *,i,itortyp(i)
1834        itortyp(i)=itype2loc(i)
1835       enddo
1836       write (iout,*) &
1837       "Val-tor parameters calculated from cumulant coefficients ntortyp"&
1838       ,ntortyp
1839       do i=-ntortyp+1,ntortyp-1
1840         do j=-ntortyp+1,ntortyp-1
1841           nterm_kcc(j,i)=2
1842           nterm_kcc_Tb(j,i)=3
1843           do k=1,nterm_kcc_Tb(j,i)
1844             do l=1,nterm_kcc_Tb(j,i)
1845               v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
1846                               +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1847               v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
1848                               +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1849             enddo
1850           enddo
1851           do k=1,nterm_kcc_Tb(j,i)
1852             do l=1,nterm_kcc_Tb(j,i)
1853 #ifdef CORRCD
1854               v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
1855                               -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1856               v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
1857                               +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1858 #else
1859               v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
1860                               -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1861               v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
1862                               +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1863 #endif
1864             enddo
1865           enddo
1866 !c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
1867         enddo
1868       enddo
1869 #else
1870       write (iout,*) "TOR_MODE>1 only with NEWCORR"
1871       stop
1872 #endif
1873       ENDIF ! TOR_MODE
1874       if (tor_mode.gt.0 .and. lprint) then
1875 !c Print valence-torsional parameters
1876         write (iout,'(a)') &
1877          "Parameters of the valence-torsional potentials"
1878         do i=-ntortyp+1,ntortyp-1
1879         do j=-ntortyp+1,ntortyp-1
1880         write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
1881         write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
1882         do k=1,nterm_kcc(j,i)
1883           do l=1,nterm_kcc_Tb(j,i)
1884             do ll=1,nterm_kcc_Tb(j,i)
1885                write (iout,'(3i5,2f15.4)')&
1886                  k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1887             enddo
1888           enddo
1889         enddo
1890         enddo
1891         enddo
1892       endif
1893
1894 #endif
1895 !elwrite(iout,*) "parmread kontrol sc-bb"
1896 ! Read of Side-chain backbone correlation parameters
1897 ! Modified 11 May 2012 by Adasko
1898 !CC
1899 !
1900      read (isccor,*) nsccortyp
1901
1902      maxinter=3
1903 !c maxinter is maximum interaction sites
1904 !write(iout,*)"maxterm_sccor",maxterm_sccor
1905 !el from module energy-------------
1906       allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1907       allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1908       allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1909       allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp))   !(maxterm_sccor,20,20)
1910 !-----------------------------------
1911       allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1912 !-----------------------------------
1913       allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1914       allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1915                -nsccortyp:nsccortyp))
1916       allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1917                -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1918       allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1919                -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1920 !-----------------------------------
1921       do i=-nsccortyp,nsccortyp
1922         do j=-nsccortyp,nsccortyp
1923           nterm_sccor(j,i)=0
1924         enddo
1925       enddo
1926 !-----------------------------------
1927
1928       read (isccor,*) (isccortyp(i),i=1,ntyp)
1929       do i=-ntyp,-1
1930         isccortyp(i)=-isccortyp(-i)
1931       enddo
1932       iscprol=isccortyp(20)
1933 !      write (iout,*) 'ntortyp',ntortyp
1934 !      maxinter=3
1935 !c maxinter is maximum interaction sites
1936       do l=1,maxinter
1937       do i=1,nsccortyp
1938         do j=1,nsccortyp
1939           read (isccor,*) &
1940       nterm_sccor(i,j),nlor_sccor(i,j)
1941           v0ijsccor=0.0d0
1942           v0ijsccor1=0.0d0
1943           v0ijsccor2=0.0d0
1944           v0ijsccor3=0.0d0
1945           si=-1.0d0
1946           nterm_sccor(-i,j)=nterm_sccor(i,j)
1947           nterm_sccor(-i,-j)=nterm_sccor(i,j)
1948           nterm_sccor(i,-j)=nterm_sccor(i,j)
1949           do k=1,nterm_sccor(i,j)
1950             read (isccor,*) kk,v1sccor(k,l,i,j),&
1951             v2sccor(k,l,i,j)
1952             if (j.eq.iscprol) then
1953              if (i.eq.isccortyp(10)) then
1954              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1955              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1956              else
1957              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1958                               +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1959              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1960                               +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1961              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1962              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1963              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1964              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1965              endif
1966             else
1967              if (i.eq.isccortyp(10)) then
1968              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1969              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1970              else
1971                if (j.eq.isccortyp(10)) then
1972              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1973              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1974                else
1975              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1976              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1977              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1978              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1979              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1980              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1981                 endif
1982                endif
1983             endif
1984             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1985             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1986             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1987             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1988             si=-si
1989            enddo
1990           do k=1,nlor_sccor(i,j)
1991             read (isccor,*) kk,vlor1sccor(k,i,j),&
1992               vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1993             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1994       (1+vlor3sccor(k,i,j)**2)
1995           enddo
1996           v0sccor(l,i,j)=v0ijsccor
1997           v0sccor(l,-i,j)=v0ijsccor1
1998           v0sccor(l,i,-j)=v0ijsccor2
1999           v0sccor(l,-i,-j)=v0ijsccor3
2000           enddo
2001         enddo
2002       enddo
2003       close (isccor)
2004       if (lprint) then
2005         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
2006         do i=1,nsccortyp
2007           do j=1,nsccortyp
2008             write (iout,*) 'ityp',i,' jtyp',j
2009             write (iout,*) 'Fourier constants'
2010             do k=1,nterm_sccor(i,j)
2011               write (iout,'(2(1pe15.5))') &
2012          (v1sccor(k,l,i,j),v2sccor(k,l,i,j),l=1,maxinter)
2013             enddo
2014             write (iout,*) 'Lorenz constants'
2015             do k=1,nlor_sccor(i,j)
2016               write (iout,'(3(1pe15.5))') &
2017                vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2018             enddo
2019           enddo
2020         enddo
2021       endif
2022
2023 !
2024 ! Read electrostatic-interaction parameters
2025 !
2026       if (lprint) then
2027         write (iout,'(/a)') 'Electrostatic interaction constants:'
2028         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2029                   'IT','JT','APP','BPP','AEL6','AEL3'
2030       endif
2031       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
2032       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
2033       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
2034       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
2035       close (ielep)
2036       do i=1,2
2037         do j=1,2
2038         rri=rpp(i,j)**6
2039         app (i,j)=epp(i,j)*rri*rri
2040         bpp (i,j)=-2.0D0*epp(i,j)*rri
2041         ael6(i,j)=elpp6(i,j)*4.2D0**6
2042         ael3(i,j)=elpp3(i,j)*4.2D0**3
2043         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2044                           ael6(i,j),ael3(i,j)
2045         enddo
2046       enddo
2047 !
2048 ! Read side-chain interaction parameters.
2049 !
2050 !el from module energy - COMMON.INTERACT-------
2051       allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2052       allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2053       allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2054       allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2055       allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2056       allocate(epslip(ntyp,ntyp))
2057       do i=1,ntyp
2058         do j=1,ntyp
2059           augm(i,j)=0.0D0
2060         enddo
2061         chip(i)=0.0D0
2062         alp(i)=0.0D0
2063         sigma0(i)=0.0D0
2064         sigii(i)=0.0D0
2065         rr0(i)=0.0D0
2066       enddo
2067 !--------------------------------
2068
2069       read (isidep,*) ipot,expon
2070 !el      if (ipot.lt.1 .or. ipot.gt.5) then
2071 !        write (iout,'(2a)') 'Error while reading SC interaction',&
2072 !                     'potential file - unknown potential type.'
2073 !        stop
2074 !wl      endif
2075       expon2=expon/2
2076       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2077        ', exponents are ',expon,2*expon
2078 !      goto (10,20,30,30,40) ipot
2079       select case(ipot)
2080 !----------------------- LJ potential ---------------------------------
2081        case (1)
2082 !   10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
2083         read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
2084         if (lprint) then
2085           write (iout,'(/a/)') 'Parameters of the LJ potential:'
2086           write (iout,'(a/)') 'The epsilon array:'
2087           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2088           write (iout,'(/a)') 'One-body parameters:'
2089           write (iout,'(a,4x,a)') 'residue','sigma'
2090           write (iout,'(a3,6x,f10.5)') (restyp(i,molnum(i)),sigma0(i),i=1,ntyp)
2091         endif
2092 !      goto 50
2093 !----------------------- LJK potential --------------------------------
2094        case (2)
2095 !   20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2096         read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2097           (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2098         if (lprint) then
2099           write (iout,'(/a/)') 'Parameters of the LJK potential:'
2100           write (iout,'(a/)') 'The epsilon array:'
2101           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2102           write (iout,'(/a)') 'One-body parameters:'
2103           write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
2104           write (iout,'(a3,6x,2f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),&
2105                 i=1,ntyp)
2106         endif
2107 !      goto 50
2108 !---------------------- GB or BP potential -----------------------------
2109        case (3:4)
2110 !   30 do i=1,ntyp
2111         if (scelemode.eq.0) then
2112         do i=1,ntyp
2113          read (isidep,*)(eps(i,j),j=i,ntyp)
2114         enddo
2115         read (isidep,*)(sigma0(i),i=1,ntyp)
2116         read (isidep,*)(sigii(i),i=1,ntyp)
2117         read (isidep,*)(chip(i),i=1,ntyp)
2118         read (isidep,*)(alp(i),i=1,ntyp)
2119         do i=1,ntyp
2120          read (isidep,*)(epslip(i,j),j=i,ntyp)
2121         enddo
2122 ! For the GB potential convert sigma'**2 into chi'
2123         if (ipot.eq.4) then
2124           do i=1,ntyp
2125             chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2126           enddo
2127         endif
2128         if (lprint) then
2129           write (iout,'(/a/)') 'Parameters of the BP potential:'
2130           write (iout,'(a/)') 'The epsilon array:'
2131           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2132           write (iout,'(/a)') 'One-body parameters:'
2133           write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',&
2134                '    chip  ','    alph  '
2135           write (iout,'(a3,6x,4f10.5)') (restyp(i,molnum(i)),sigma0(i),sigii(i),&
2136                            chip(i),alp(i),i=1,ntyp)
2137         endif
2138         else
2139       allocate(icharge(ntyp1))
2140 !      print *,ntyp,icharge(i)
2141       icharge(:)=0
2142       read (isidep,*) (icharge(i),i=1,ntyp)
2143       print *,ntyp,icharge(i)
2144 !      if(.not.allocated(eps)) allocate(eps(-ntyp
2145       write (2,*) "icharge",(icharge(i),i=1,ntyp)
2146        allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2147        allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2148        allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2149        allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2150        allocate(epsintab(ntyp,ntyp))
2151        allocate(dtail(2,ntyp,ntyp))
2152       print *,"control line 1"
2153        allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2154        allocate(wqdip(2,ntyp,ntyp))
2155        allocate(wstate(4,ntyp,ntyp))
2156        allocate(dhead(2,2,ntyp,ntyp))
2157        allocate(nstate(ntyp,ntyp))
2158        allocate(debaykap(ntyp,ntyp))
2159       print *,"control line 2"
2160       if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2161       if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2162
2163       do i=1,ntyp
2164        do j=1,i
2165 !        write (*,*) "Im in ALAB", i, " ", j
2166         read(isidep,*) &
2167        eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), &
2168        (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &
2169        chis(i,j),chis(j,i), &
2170        nstate(i,j),(wstate(k,i,j),k=1,4), &
2171        dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),&
2172        dtail(1,i,j),dtail(2,i,j), &
2173        epshead(i,j),sig0head(i,j), &
2174        rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), &
2175        alphapol(i,j),alphapol(j,i), &
2176        (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j)
2177 !       print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
2178        END DO
2179       END DO
2180       DO i = 1, ntyp
2181        DO j = i+1, ntyp
2182         eps(i,j) = eps(j,i)
2183         sigma(i,j) = sigma(j,i)
2184         sigmap1(i,j)=sigmap1(j,i)
2185         sigmap2(i,j)=sigmap2(j,i)
2186         sigiso1(i,j)=sigiso1(j,i)
2187         sigiso2(i,j)=sigiso2(j,i)
2188 !        print *,"ATU",sigma(j,i),sigma(i,j),i,j
2189         nstate(i,j) = nstate(j,i)
2190         dtail(1,i,j) = dtail(1,j,i)
2191         dtail(2,i,j) = dtail(2,j,i)
2192         DO k = 1, 4
2193          alphasur(k,i,j) = alphasur(k,j,i)
2194          wstate(k,i,j) = wstate(k,j,i)
2195          alphiso(k,i,j) = alphiso(k,j,i)
2196         END DO
2197
2198         dhead(2,1,i,j) = dhead(1,1,j,i)
2199         dhead(2,2,i,j) = dhead(1,2,j,i)
2200         dhead(1,1,i,j) = dhead(2,1,j,i)
2201         dhead(1,2,i,j) = dhead(2,2,j,i)
2202
2203         epshead(i,j) = epshead(j,i)
2204         sig0head(i,j) = sig0head(j,i)
2205
2206         DO k = 1, 2
2207          wqdip(k,i,j) = wqdip(k,j,i)
2208         END DO
2209
2210         wquad(i,j) = wquad(j,i)
2211         epsintab(i,j) = epsintab(j,i)
2212         debaykap(i,j)=debaykap(j,i)
2213 !        if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2214        END DO
2215       END DO
2216       endif
2217
2218 !      goto 50
2219 !--------------------- GBV potential -----------------------------------
2220        case (5)
2221 !   40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2222         read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2223           (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2224         (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2225         if (lprint) then
2226           write (iout,'(/a/)') 'Parameters of the GBV potential:'
2227           write (iout,'(a/)') 'The epsilon array:'
2228           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2229           write (iout,'(/a)') 'One-body parameters:'
2230           write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',&
2231             's||/s_|_^2','    chip  ','    alph  '
2232           write (iout,'(a3,6x,5f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),&
2233                  sigii(i),chip(i),alp(i),i=1,ntyp)
2234         endif
2235        case default
2236         write (iout,'(2a)') 'Error while reading SC interaction',&
2237                      'potential file - unknown potential type.'
2238         stop
2239 !   50 continue
2240       end select
2241 !      continue
2242       close (isidep)
2243 !-----------------------------------------------------------------------
2244 ! Calculate the "working" parameters of SC interactions.
2245
2246 !el from module energy - COMMON.INTERACT-------
2247 !      allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2248             if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1))
2249       allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1)) !(ntyp,ntyp)
2250       allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2251       if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2252       allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2253       allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2254         dcavtub(ntyp1))
2255       allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2256         tubetranene(ntyp1))
2257       do i=1,ntyp1
2258         do j=1,ntyp1
2259           aa_aq(i,j)=0.0D0
2260           bb_aq(i,j)=0.0D0
2261           aa_lip(i,j)=0.0d0
2262           bb_lip(i,j)=0.0d0
2263                if (scelemode.eq.0) then
2264           chi(i,j)=0.0D0
2265           sigma(i,j)=0.0D0
2266           r0(i,j)=0.0D0
2267             endif
2268         enddo
2269       enddo
2270 !--------------------------------
2271
2272       do i=2,ntyp
2273         do j=1,i-1
2274           eps(i,j)=eps(j,i)
2275           epslip(i,j)=epslip(j,i)
2276         enddo
2277       enddo
2278       if (scelemode.eq.0) then
2279       do i=1,ntyp
2280         do j=i,ntyp
2281           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2282           sigma(j,i)=sigma(i,j)
2283           rs0(i,j)=dwa16*sigma(i,j)
2284           rs0(j,i)=rs0(i,j)
2285         enddo
2286       enddo
2287       endif
2288       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2289        'Working parameters of the SC interactions:',&
2290        '     a    ','     b    ','   augm   ','  sigma ','   r0   ',&
2291        '  chi1   ','   chi2   '
2292       do i=1,ntyp
2293         do j=i,ntyp
2294           epsij=eps(i,j)
2295           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2296             rrij=sigma(i,j)
2297             print *,"rrij",rrij
2298           else
2299             rrij=rr0(i)+rr0(j)
2300           endif
2301           r0(i,j)=rrij
2302           r0(j,i)=rrij
2303           rrij=rrij**expon
2304           epsij=eps(i,j)
2305           sigeps=dsign(1.0D0,epsij)
2306           epsij=dabs(epsij)
2307           aa_aq(i,j)=epsij*rrij*rrij
2308           bb_aq(i,j)=-sigeps*epsij*rrij
2309           aa_aq(j,i)=aa_aq(i,j)
2310           bb_aq(j,i)=bb_aq(i,j)
2311           epsijlip=epslip(i,j)
2312           sigeps=dsign(1.0D0,epsijlip)
2313           epsijlip=dabs(epsijlip)
2314           aa_lip(i,j)=epsijlip*rrij*rrij
2315           bb_lip(i,j)=-sigeps*epsijlip*rrij
2316           aa_lip(j,i)=aa_lip(i,j)
2317           bb_lip(j,i)=bb_lip(i,j)
2318           if ((ipot.gt.2).and. (scelemode.eq.0))then
2319             sigt1sq=sigma0(i)**2
2320             sigt2sq=sigma0(j)**2
2321             sigii1=sigii(i)
2322             sigii2=sigii(j)
2323             ratsig1=sigt2sq/sigt1sq
2324             ratsig2=1.0D0/ratsig1
2325             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2326             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2327             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2328           else
2329             rsum_max=sigma(i,j)
2330           endif
2331 !         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2332             sigmaii(i,j)=rsum_max
2333             sigmaii(j,i)=rsum_max
2334 !         else
2335 !           sigmaii(i,j)=r0(i,j)
2336 !           sigmaii(j,i)=r0(i,j)
2337 !         endif
2338 !d        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2339           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2340             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2341             augm(i,j)=epsij*r_augm**(2*expon)
2342 !           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2343             augm(j,i)=augm(i,j)
2344           else
2345             augm(i,j)=0.0D0
2346             augm(j,i)=0.0D0
2347           endif
2348           if (lprint) then
2349             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')  &
2350             restyp(i,molnum(i)),restyp(j,molnum(i)),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2351             sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2352           endif
2353         enddo
2354       enddo
2355       allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2356       allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2357       allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2358       allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2359       allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2360       allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2361       allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2362       allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2363       allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2364       allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2365       allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2366       allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2367       allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2368       allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2369       allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2370       allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2371
2372 !      augm(:,:)=0.0D0
2373 !      chip(:)=0.0D0
2374 !      alp(:)=0.0D0
2375 !      sigma0(:)=0.0D0
2376 !      sigii(:)=0.0D0
2377 !      rr0(:)=0.0D0
2378
2379       read (isidep_nucl,*) ipot_nucl
2380 !      print *,"TU?!",ipot_nucl
2381       if (ipot_nucl.eq.1) then
2382         do i=1,ntyp_molec(2)
2383           do j=i,ntyp_molec(2)
2384             read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2385             elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2386           enddo
2387         enddo
2388       else
2389         do i=1,ntyp_molec(2)
2390           do j=i,ntyp_molec(2)
2391             read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2392                chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2393                elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2394           enddo
2395         enddo
2396       endif
2397 !      rpp(1,1)=2**(1.0/6.0)*5.16158
2398       do i=1,ntyp_molec(2)
2399         do j=i,ntyp_molec(2)
2400           rrij=sigma_nucl(i,j)
2401           r0_nucl(i,j)=rrij
2402           r0_nucl(j,i)=rrij
2403           rrij=rrij**expon
2404           epsij=4*eps_nucl(i,j)
2405           sigeps=dsign(1.0D0,epsij)
2406           epsij=dabs(epsij)
2407           aa_nucl(i,j)=epsij*rrij*rrij
2408           bb_nucl(i,j)=-sigeps*epsij*rrij
2409           ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2410           ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2411           ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2412           ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2413           sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2414          2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2415         enddo
2416         do j=1,i-1
2417           aa_nucl(i,j)=aa_nucl(j,i)
2418           bb_nucl(i,j)=bb_nucl(j,i)
2419           ael3_nucl(i,j)=ael3_nucl(j,i)
2420           ael6_nucl(i,j)=ael6_nucl(j,i)
2421           ael63_nucl(i,j)=ael63_nucl(j,i)
2422           ael32_nucl(i,j)=ael32_nucl(j,i)
2423           elpp3_nucl(i,j)=elpp3_nucl(j,i)
2424           elpp6_nucl(i,j)=elpp6_nucl(j,i)
2425           elpp63_nucl(i,j)=elpp63_nucl(j,i)
2426           elpp32_nucl(i,j)=elpp32_nucl(j,i)
2427           eps_nucl(i,j)=eps_nucl(j,i)
2428           sigma_nucl(i,j)=sigma_nucl(j,i)
2429           sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2430         enddo
2431       enddo
2432
2433       write(iout,*) "tube param"
2434       read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2435       ccavtubpep,dcavtubpep,tubetranenepep
2436       sigmapeptube=sigmapeptube**6
2437       sigeps=dsign(1.0D0,epspeptube)
2438       epspeptube=dabs(epspeptube)
2439       pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2440       pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2441       write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2442       do i=1,ntyp
2443        read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2444       ccavtub(i),dcavtub(i),tubetranene(i)
2445        sigmasctube=sigmasctube**6
2446        sigeps=dsign(1.0D0,epssctube)
2447        epssctube=dabs(epssctube)
2448        sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2449        sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2450       write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2451       enddo
2452 !-----------------READING SC BASE POTENTIALS-----------------------------
2453       allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2454       allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2455       allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2456       allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2457       allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2458       allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2459       allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2460       allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2461       allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2462       allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2463       allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2464       allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2465       allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2466       allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2467       allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2468       allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2469
2470
2471       do i=1,ntyp_molec(1)
2472        do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2473         write (*,*) "Im in ", i, " ", j
2474         read(isidep_scbase,*) &
2475         eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2476         chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2477          write(*,*) "eps",eps_scbase(i,j)
2478         read(isidep_scbase,*) &
2479        (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2480        chis_scbase(i,j,1),chis_scbase(i,j,2)
2481         read(isidep_scbase,*) &
2482        dhead_scbasei(i,j), &
2483        dhead_scbasej(i,j), &
2484        rborn_scbasei(i,j),rborn_scbasej(i,j)
2485         read(isidep_scbase,*) &
2486        (wdipdip_scbase(k,i,j),k=1,3), &
2487        (wqdip_scbase(k,i,j),k=1,2)
2488         read(isidep_scbase,*) &
2489        alphapol_scbase(i,j), &
2490        epsintab_scbase(i,j)
2491        END DO
2492       END DO
2493       allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2494       allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2495
2496       do i=1,ntyp_molec(1)
2497        do j=1,ntyp_molec(2)-1
2498           epsij=eps_scbase(i,j)
2499           rrij=sigma_scbase(i,j)
2500 !          r0(i,j)=rrij
2501 !          r0(j,i)=rrij
2502           rrij=rrij**expon
2503 !          epsij=eps(i,j)
2504           sigeps=dsign(1.0D0,epsij)
2505           epsij=dabs(epsij)
2506           aa_scbase(i,j)=epsij*rrij*rrij
2507           bb_scbase(i,j)=-sigeps*epsij*rrij
2508         enddo
2509        enddo
2510 !-----------------READING PEP BASE POTENTIALS-------------------
2511       allocate(eps_pepbase(ntyp_molec(2)))
2512       allocate(sigma_pepbase(ntyp_molec(2)))
2513       allocate(chi_pepbase(ntyp_molec(2),2))
2514       allocate(chipp_pepbase(ntyp_molec(2),2))
2515       allocate(alphasur_pepbase(4,ntyp_molec(2)))
2516       allocate(sigmap1_pepbase(ntyp_molec(2)))
2517       allocate(sigmap2_pepbase(ntyp_molec(2)))
2518       allocate(chis_pepbase(ntyp_molec(2),2))
2519       allocate(wdipdip_pepbase(3,ntyp_molec(2)))
2520
2521
2522        do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2523         write (*,*) "Im in ", i, " ", j
2524         read(isidep_pepbase,*) &
2525         eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
2526         chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
2527          write(*,*) "eps",eps_pepbase(j)
2528         read(isidep_pepbase,*) &
2529        (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
2530        chis_pepbase(j,1),chis_pepbase(j,2)
2531         read(isidep_pepbase,*) &
2532        (wdipdip_pepbase(k,j),k=1,3)
2533        END DO
2534       allocate(aa_pepbase(ntyp_molec(2)))
2535       allocate(bb_pepbase(ntyp_molec(2)))
2536
2537        do j=1,ntyp_molec(2)-1
2538           epsij=eps_pepbase(j)
2539           rrij=sigma_pepbase(j)
2540 !          r0(i,j)=rrij
2541 !          r0(j,i)=rrij
2542           rrij=rrij**expon
2543 !          epsij=eps(i,j)
2544           sigeps=dsign(1.0D0,epsij)
2545           epsij=dabs(epsij)
2546           aa_pepbase(j)=epsij*rrij*rrij
2547           bb_pepbase(j)=-sigeps*epsij*rrij
2548         enddo
2549 !--------------READING SC PHOSPHATE-------------------------------------
2550 !--------------READING SC PHOSPHATE-------------------------------------
2551       allocate(eps_scpho(ntyp_molec(1)))
2552       allocate(sigma_scpho(ntyp_molec(1)))
2553       allocate(chi_scpho(ntyp_molec(1),2))
2554       allocate(chipp_scpho(ntyp_molec(1),2))
2555       allocate(alphasur_scpho(4,ntyp_molec(1)))
2556       allocate(sigmap1_scpho(ntyp_molec(1)))
2557       allocate(sigmap2_scpho(ntyp_molec(1)))
2558       allocate(chis_scpho(ntyp_molec(1),2))
2559       allocate(wqq_scpho(ntyp_molec(1)))
2560       allocate(wqdip_scpho(2,ntyp_molec(1)))
2561       allocate(alphapol_scpho(ntyp_molec(1)))
2562       allocate(epsintab_scpho(ntyp_molec(1)))
2563       allocate(dhead_scphoi(ntyp_molec(1)))
2564       allocate(rborn_scphoi(ntyp_molec(1)))
2565       allocate(rborn_scphoj(ntyp_molec(1)))
2566       allocate(alphi_scpho(ntyp_molec(1)))
2567
2568
2569 !      j=1
2570        do j=1,ntyp_molec(1) ! without U then we will take T for U
2571         write (*,*) "Im in scpho ", i, " ", j
2572         read(isidep_scpho,*) &
2573         eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
2574         chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
2575          write(*,*) "eps",eps_scpho(j)
2576         read(isidep_scpho,*) &
2577        (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
2578        chis_scpho(j,1),chis_scpho(j,2)
2579         read(isidep_scpho,*) &
2580        (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
2581         read(isidep_scpho,*) &
2582          epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
2583          alphi_scpho(j)
2584
2585        END DO
2586       allocate(aa_scpho(ntyp_molec(1)))
2587       allocate(bb_scpho(ntyp_molec(1)))
2588        do j=1,ntyp_molec(1)
2589           epsij=eps_scpho(j)
2590           rrij=sigma_scpho(j)
2591 !          r0(i,j)=rrij
2592 !          r0(j,i)=rrij
2593           rrij=rrij**expon
2594 !          epsij=eps(i,j)
2595           sigeps=dsign(1.0D0,epsij)
2596           epsij=dabs(epsij)
2597           aa_scpho(j)=epsij*rrij*rrij
2598           bb_scpho(j)=-sigeps*epsij*rrij
2599         enddo
2600
2601
2602         read(isidep_peppho,*) &
2603         eps_peppho,sigma_peppho
2604         read(isidep_peppho,*) &
2605        (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
2606         read(isidep_peppho,*) &
2607        (wqdip_peppho(k),k=1,2)
2608
2609           epsij=eps_peppho
2610           rrij=sigma_peppho
2611 !          r0(i,j)=rrij
2612 !          r0(j,i)=rrij
2613           rrij=rrij**expon
2614 !          epsij=eps(i,j)
2615           sigeps=dsign(1.0D0,epsij)
2616           epsij=dabs(epsij)
2617           aa_peppho=epsij*rrij*rrij
2618           bb_peppho=-sigeps*epsij*rrij
2619
2620
2621       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2622       do i=1,ntyp
2623         do j=1,2
2624           bad(i,j)=0.0D0
2625         enddo
2626       enddo
2627 #ifdef CLUSTER
2628 !
2629 ! Define the SC-p interaction constants
2630 !
2631       do i=1,20
2632         do j=1,2
2633           eps_scp(i,j)=-1.5d0
2634           rscp(i,j)=4.0d0
2635         enddo
2636       enddo
2637 #endif
2638       allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2639       read (itorp_nucl,*) ntortyp_nucl
2640 !      print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2641 !el from energy module---------
2642       allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2643       allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2644
2645       allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2646       allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2647       allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2648       allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2649
2650       allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2651       allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2652 !el---------------------------
2653       nterm_nucl(:,:)=0
2654       nlor_nucl(:,:)=0
2655 !el--------------------
2656       read (itorp_nucl,*) &
2657         (itortyp_nucl(i),i=1,ntyp_molec(2))
2658 !        print *,itortyp_nucl(:)
2659 !c      write (iout,*) 'ntortyp',ntortyp
2660       do i=1,ntortyp_nucl
2661         do j=1,ntortyp_nucl
2662           read (itorp_nucl,*) nterm_nucl(i,j),nlor_nucl(i,j)
2663 !           print *,nterm_nucl(i,j),nlor_nucl(i,j)
2664           v0ij=0.0d0
2665           si=-1.0d0
2666           do k=1,nterm_nucl(i,j)
2667             read (itorp_nucl,*) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2668             v0ij=v0ij+si*v1_nucl(k,i,j)
2669             si=-si
2670           enddo
2671           do k=1,nlor_nucl(i,j)
2672             read (itorp_nucl,*) kk,vlor1_nucl(k,i,j),&
2673               vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2674             v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2675           enddo
2676           v0_nucl(i,j)=v0ij
2677         enddo
2678       enddo
2679
2680
2681 !elwrite(iout,*) "parmread kontrol before oldscp"
2682 !
2683 ! Define the SC-p interaction constants
2684 !
2685 #ifdef OLDSCP
2686       do i=1,20
2687 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2688 ! helix formation)
2689 !       aad(i,1)=0.3D0*4.0D0**12
2690 ! Following line for constants currently implemented
2691 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2692         aad(i,1)=1.5D0*4.0D0**12
2693 !       aad(i,1)=0.17D0*5.6D0**12
2694         aad(i,2)=aad(i,1)
2695 ! "Soft" SC-p repulsion
2696         bad(i,1)=0.0D0
2697 ! Following line for constants currently implemented
2698 !       aad(i,1)=0.3D0*4.0D0**6
2699 ! "Hard" SC-p repulsion
2700         bad(i,1)=3.0D0*4.0D0**6
2701 !       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2702         bad(i,2)=bad(i,1)
2703 !       aad(i,1)=0.0D0
2704 !       aad(i,2)=0.0D0
2705 !       bad(i,1)=1228.8D0
2706 !       bad(i,2)=1228.8D0
2707       enddo
2708 #else
2709 !
2710 ! 8/9/01 Read the SC-p interaction constants from file
2711 !
2712       do i=1,ntyp
2713         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
2714       enddo
2715       do i=1,ntyp
2716         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2717         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2718         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2719         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2720       enddo
2721
2722       if (lprint) then
2723         write (iout,*) "Parameters of SC-p interactions:"
2724         do i=1,20
2725           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2726            eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2727         enddo
2728       endif
2729 #endif
2730       allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
2731
2732       do i=1,ntyp_molec(2)
2733         read (iscpp_nucl,*) eps_scp_nucl(i),rscp_nucl(i)
2734       enddo
2735       do i=1,ntyp_molec(2)
2736         aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
2737         bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
2738       enddo
2739       r0pp=1.12246204830937298142*5.16158
2740       epspp=4.95713/4
2741       AEES=108.661
2742       BEES=0.433246
2743
2744 !
2745 ! Define the constants of the disulfide bridge
2746 !
2747       ebr=-5.50D0
2748 !
2749 ! Old arbitrary potential - commented out.
2750 !
2751 !      dbr= 4.20D0
2752 !      fbr= 3.30D0
2753 !
2754 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2755 ! energy surface of diethyl disulfide.
2756 ! A. Liwo and U. Kozlowska, 11/24/03
2757 !
2758       D0CM = 3.78d0
2759       AKCM = 15.1d0
2760       AKTH = 11.0d0
2761       AKCT = 12.0d0
2762       V1SS =-1.08d0
2763       V2SS = 7.61d0
2764       V3SS = 13.7d0
2765
2766       if (lprint) then
2767       write (iout,'(/a)') "Disulfide bridge parameters:"
2768       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2769       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2770       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2771       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2772        ' v3ss:',v3ss
2773       endif
2774       return
2775       end subroutine parmread
2776 #ifndef CLUSTER
2777 !-----------------------------------------------------------------------------
2778 ! mygetenv.F
2779 !-----------------------------------------------------------------------------
2780       subroutine mygetenv(string,var)
2781 !
2782 ! Version 1.0
2783 !
2784 ! This subroutine passes the environmental variables to FORTRAN program.
2785 ! If the flags -DMYGETENV and -DMPI are not for compilation, it calls the
2786 ! standard FORTRAN GETENV subroutine. If both flags are set, the subroutine
2787 ! reads the environmental variables from $HOME/.env
2788 !
2789 ! Usage: As for the standard FORTRAN GETENV subroutine.
2790 !
2791 ! Purpose: some versions/installations of MPI do not transfer the environmental
2792 ! variables to slave processors, if these variables are set in the shell script
2793 ! from which mpirun is called.
2794 !
2795 ! A.Liwo, 7/29/01
2796 !
2797 #ifdef MPI
2798       use MPI_data
2799       include "mpif.h"
2800 #endif
2801 !      implicit none
2802       character*(*) :: string,var
2803 #if defined(MYGETENV) && defined(MPI)
2804 !      include "DIMENSIONS.ZSCOPT"
2805 !      include "mpif.h"
2806 !      include "COMMON.MPI"
2807 !el      character*360 ucase
2808 !el      external ucase
2809       character(len=360) :: string1(360),karta
2810       character(len=240) :: home
2811       integer i,n !,ilen
2812 !el      external ilen
2813       call getenv("HOME",home)
2814       open(99,file=home(:ilen(home))//"/.env",status="OLD",err=112)
2815       do while (.true.)
2816         read (99,end=111,err=111,'(a)') karta
2817         do i=1,80
2818           string1(i)=" "
2819         enddo
2820         call split_string(karta,string1,80,n)
2821         if (ucase(string1(1)(:ilen(string1(1)))).eq."SETENV" .and. &
2822          string1(2)(:ilen(string1(2))).eq.string(:ilen(string)) ) then
2823            var=string1(3)
2824            print *,"Processor",me,": ",var(:ilen(var)),&
2825             " assigned to ",string(:ilen(string))
2826            close(99)
2827            return
2828         endif
2829       enddo
2830  111  print *,"Environment variable ",string(:ilen(string))," not set."
2831       close(99)
2832       return
2833  112  print *,"Error opening environment file!"
2834 #else
2835       call getenv(string,var)
2836 #endif
2837       return
2838       end subroutine mygetenv
2839 !-----------------------------------------------------------------------------
2840 ! readrtns.F
2841 !-----------------------------------------------------------------------------
2842       subroutine read_general_data(*)
2843
2844       use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions,&
2845           scelemode,TUBEmode,tor_mode
2846
2847       use energy_data, only:distchainmax,tubeR0,tubecenter
2848       use geometry_data, only:boxxsize,boxysize,boxzsize,bordtubetop,&
2849           bordtubebot,tubebufthick,buftubebot,buftubetop
2850 !      implicit none
2851 !      include "DIMENSIONS"
2852 !      include "DIMENSIONS.ZSCOPT"
2853 !      include "DIMENSIONS.FREE"
2854 !      include "COMMON.TORSION"
2855 !      include "COMMON.INTERACT"
2856 !      include "COMMON.IOUNITS"
2857 !      include "COMMON.TIME1"
2858 !      include "COMMON.PROT"
2859 !      include "COMMON.PROTFILES"
2860 !      include "COMMON.CHAIN"
2861 !      include "COMMON.NAMES"
2862 !      include "COMMON.FFIELD"
2863 !      include "COMMON.ENEPS"
2864 !      include "COMMON.WEIGHTS"
2865 !      include "COMMON.FREE"
2866 !      include "COMMON.CONTROL"
2867 !      include "COMMON.ENERGIES"
2868       character(len=800) :: controlcard
2869       integer :: i,j,k,ii,n_ene_found
2870       integer :: ind,itype1,itype2,itypf,itypsc,itypp
2871 !el      integer ilen
2872 !el      external ilen
2873 !el      character*16 ucase
2874       character(len=16) :: key
2875 !el      external ucase
2876       call card_concat(controlcard,.true.)
2877       call readi(controlcard,"N_ENE",n_eneW,max_eneW)
2878       if (n_eneW.gt.max_eneW) then
2879         write (iout,*) "Error: parameter out of range: N_ENE",n_eneW,&
2880           max_eneW
2881         return 1
2882       endif
2883       call readi(controlcard,"NPARMSET",nparmset,1)
2884 !elwrite(iout,*)"in read_gen data"
2885       separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
2886       call readi(controlcard,"IPARMPRINT",iparmprint,1)
2887       write (iout,*) "PARMPRINT",iparmprint
2888       if (nparmset.gt.max_parm) then
2889         write (iout,*) "Error: parameter out of range: NPARMSET",&
2890           nparmset, Max_Parm
2891         return 1
2892       endif
2893 !elwrite(iout,*)"in read_gen data"
2894       call readi(controlcard,"MAXIT",maxit,5000)
2895       call reada(controlcard,"FIMIN",fimin,1.0d-3)
2896       call readi(controlcard,"ENSEMBLES",ensembles,0)
2897       hamil_rep=index(controlcard,"HAMIL_REP").gt.0
2898       write (iout,*) "Number of energy parameter sets",nparmset
2899       allocate(isampl(nparmset))
2900       call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
2901       write (iout,*) "MaxSlice",MaxSlice
2902       call readi(controlcard,"NSLICE",nslice,1)
2903 !elwrite(iout,*)"in read_gen data"
2904       call flush(iout)
2905       if (nslice.gt.MaxSlice) then
2906         write (iout,*) "Error: parameter out of range: NSLICE",nslice,&
2907           MaxSlice
2908         return 1
2909       endif
2910       write (iout,*) "Frequency of storing conformations",&
2911        (isampl(i),i=1,nparmset)
2912       write (iout,*) "Maxit",maxit," Fimin",fimin
2913       call readi(controlcard,"NQ",nQ,1)
2914       if (nQ.gt.MaxQ) then
2915         write (iout,*) "Error: parameter out of range: NQ",nq,&
2916           maxq
2917         return 1
2918       endif
2919       indpdb=0
2920       if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
2921       call reada(controlcard,"DELTA",delta,1.0d-2)
2922       call readi(controlcard,"EINICHECK",einicheck,2)
2923       call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
2924       call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
2925       call readi(controlcard,"RESCALE",rescale_modeW,1)
2926       call reada(controlcard,'BOXX',boxxsize,100.0d0)
2927       call reada(controlcard,'BOXY',boxysize,100.0d0)
2928       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
2929       call readi(controlcard,"SCELEMODE",scelemode,0)
2930       print *,"SCELE",scelemode
2931       call readi(controlcard,'TORMODE',tor_mode,0)
2932 !C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2933         write(iout,*) "torsional and valence angle mode",tor_mode
2934
2935       call readi(controlcard,'TUBEMOD',tubemode,0)
2936
2937       if (TUBEmode.gt.0) then
2938        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
2939        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
2940        call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
2941        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
2942        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
2943        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
2944        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
2945        buftubebot=bordtubebot+tubebufthick
2946        buftubetop=bordtubetop-tubebufthick
2947       endif
2948       ions=index(controlcard,"IONS").gt.0
2949       call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
2950       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
2951       write(iout,*) "R_CUT_ELE=",r_cut_ele
2952       check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
2953       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
2954       call readi(controlcard,'SYM',symetr,1)
2955       write (iout,*) "DISTCHAINMAX",distchainmax
2956       write (iout,*) "delta",delta
2957       write (iout,*) "einicheck",einicheck
2958       write (iout,*) "rescale_mode",rescale_modeW
2959       call flush(iout)
2960       bxfile=index(controlcard,"BXFILE").gt.0
2961       cxfile=index(controlcard,"CXFILE").gt.0
2962       if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile) &
2963        bxfile=.true.
2964       histfile=index(controlcard,"HISTFILE").gt.0
2965       histout=index(controlcard,"HISTOUT").gt.0
2966       entfile=index(controlcard,"ENTFILE").gt.0
2967       zscfile=index(controlcard,"ZSCFILE").gt.0
2968       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
2969       write (iout,*) "with_dihed_constr ",with_dihed_constr
2970       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
2971       return
2972       end subroutine read_general_data
2973 !------------------------------------------------------------------------------
2974       subroutine read_efree(*)
2975 !
2976 ! Read molecular data
2977 !
2978 !      implicit none
2979 !      include 'DIMENSIONS'
2980 !      include 'DIMENSIONS.ZSCOPT'
2981 !      include 'DIMENSIONS.COMPAR'
2982 !      include 'DIMENSIONS.FREE'
2983 !      include 'COMMON.IOUNITS'
2984 !      include 'COMMON.TIME1'
2985 !      include 'COMMON.SBRIDGE'
2986 !      include 'COMMON.CONTROL'
2987 !      include 'COMMON.CHAIN'
2988 !      include 'COMMON.HEADER'
2989 !      include 'COMMON.GEO'
2990 !      include 'COMMON.FREE'
2991       character(len=320) :: controlcard !,ucase
2992       integer :: iparm,ib,i,j,npars
2993 !el      integer ilen
2994 !el      external ilen
2995
2996       if (hamil_rep) then
2997         npars=1
2998       else
2999         npars=nParmSet
3000       endif
3001
3002 !      call alloc_wham_arrays
3003 !      allocate(nT_h(nParmSet))
3004 !      allocate(replica(nParmSet))
3005 !      allocate(umbrella(nParmSet))
3006 !      allocate(read_iset(nParmSet))
3007 !      allocate(nT_h(nParmSet))
3008
3009       do iparm=1,npars
3010
3011       call card_concat(controlcard,.true.)
3012       call readi(controlcard,'NT',nT_h(iparm),1)
3013       write (iout,*) "iparm",iparm," nt",nT_h(iparm)
3014       call flush(iout)
3015       if (nT_h(iparm).gt.MaxT_h) then
3016         write (iout,*)  "Error: parameter out of range: NT",nT_h(iparm),&
3017           MaxT_h
3018         return 1
3019       endif
3020       replica(iparm)=index(controlcard,"REPLICA").gt.0
3021       umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
3022       read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
3023       write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",&
3024         replica(iparm)," umbrella ",umbrella(iparm),&
3025         " read_iset",read_iset(iparm)
3026       call flush(iout)
3027       do ib=1,nT_h(iparm)
3028         call card_concat(controlcard,.true.)
3029         call readi(controlcard,'NR',nR(ib,iparm),1)
3030         if (umbrella(iparm)) then
3031           nRR(ib,iparm)=1
3032         else
3033           nRR(ib,iparm)=nR(ib,iparm)
3034         endif
3035         if (nR(ib,iparm).gt.MaxR) then
3036           write (iout,*)  "Error: parameter out of range: NR",&
3037             nR(ib,iparm),MaxR
3038         return 1
3039         endif
3040         call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
3041         beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
3042         call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),&
3043           0.0d0)
3044         do i=1,nR(ib,iparm)
3045           call card_concat(controlcard,.true.)
3046           call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,&
3047             100.0d0)
3048           call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,&
3049             0.0d0)
3050         enddo
3051       enddo
3052       do ib=1,nT_h(iparm)
3053         write (iout,*) "ib",ib," beta_h",&
3054           1.0d0/(0.001987*beta_h(ib,iparm))
3055         write (iout,*) "nR",nR(ib,iparm)
3056         write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
3057         do i=1,nR(ib,iparm)
3058           write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),&
3059             "q0",(q0(j,i,ib,iparm),j=1,nQ)
3060         enddo
3061         call flush(iout)
3062       enddo
3063
3064       enddo
3065
3066       if (hamil_rep) then
3067
3068        do iparm=2,nParmSet
3069           nT_h(iparm)=nT_h(1)
3070          do ib=1,nT_h(iparm)
3071            nR(ib,iparm)=nR(ib,1)
3072            if (umbrella(iparm)) then
3073              nRR(ib,iparm)=1
3074            else
3075              nRR(ib,iparm)=nR(ib,1)
3076            endif
3077            beta_h(ib,iparm)=beta_h(ib,1)
3078            do i=1,nR(ib,iparm)
3079              f(i,ib,iparm)=f(i,ib,1)
3080              do j=1,nQ
3081                KH(j,i,ib,iparm)=KH(j,i,ib,1)
3082                Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
3083              enddo
3084            enddo
3085            replica(iparm)=replica(1)
3086            umbrella(iparm)=umbrella(1)
3087            read_iset(iparm)=read_iset(1)
3088          enddo
3089        enddo
3090
3091       endif
3092
3093       return
3094       end subroutine read_efree
3095 !-----------------------------------------------------------------------------
3096       subroutine read_protein_data(*)
3097 !      implicit none
3098 !      include "DIMENSIONS"
3099 !      include "DIMENSIONS.ZSCOPT"
3100 !      include "DIMENSIONS.FREE"
3101 #ifdef MPI
3102       use MPI_data
3103       include "mpif.h"
3104       integer :: IERROR,ERRCODE!,STATUS(MPI_STATUS_SIZE)
3105 !      include "COMMON.MPI"
3106 #endif
3107 !      include "COMMON.CHAIN"
3108 !      include "COMMON.IOUNITS"
3109 !      include "COMMON.PROT"
3110 !      include "COMMON.PROTFILES"
3111 !      include "COMMON.NAMES"
3112 !      include "COMMON.FREE"
3113 !      include "COMMON.OBCINKA"
3114       character(len=64) :: nazwa
3115       character(len=16000) :: controlcard
3116       integer :: i,ii,ib,iR,iparm,nthr,npars !,ilen,iroof
3117 !el      external ilen,iroof
3118       if (hamil_rep) then
3119         npars=1
3120       else
3121         npars=nparmset
3122       endif
3123
3124       do iparm=1,npars
3125
3126 ! Read names of files with conformation data.
3127       if (replica(iparm)) then
3128         nthr = 1
3129       else
3130         nthr = nT_h(iparm)
3131       endif
3132       do ib=1,nthr
3133       do ii=1,nRR(ib,iparm)
3134       write (iout,*) "Parameter set",iparm," temperature",ib,&
3135        " window",ii
3136       call flush(iout)
3137       call card_concat(controlcard,.true.)
3138       write (iout,*) controlcard(:ilen(controlcard))
3139       call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
3140       call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
3141       call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
3142       call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
3143       call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),&
3144        maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
3145       call reada(controlcard,"TIME_START",&
3146         time_start_collect(ii,ib,iparm),0.0d0)
3147       call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),&
3148         1.0d10)
3149       write (iout,*) "rec_start",rec_start(ii,ib,iparm),&
3150        " rec_end",rec_end(ii,ib,iparm)
3151       write (iout,*) "time_start",time_start_collect(ii,ib,iparm),&
3152        " time_end",time_end_collect(ii,ib,iparm)
3153       call flush(iout)
3154       if (replica(iparm)) then
3155         call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
3156         write (iout,*) "Number of trajectories",totraj(ii,iparm)
3157         call flush(iout)
3158       endif
3159       if (nfile_bin(ii,ib,iparm).lt.2 &
3160           .and. nfile_asc(ii,ib,iparm).eq.0 &
3161           .and. nfile_cx(ii,ib,iparm).eq.0) then
3162         write (iout,*) "Error - no action specified!"
3163         return 1
3164       endif
3165       if (nfile_bin(ii,ib,iparm).gt.0) then
3166         call card_concat(controlcard,.false.)
3167         call split_string(controlcard,protfiles(1,1,ii,ib,iparm),&
3168          maxfile_prot,nfile_bin(ii,ib,iparm))
3169 #ifdef DEBUG
3170         write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
3171         write(iout,*) (protfiles(i,1,ii,ib,iparm),&
3172           i=1,nfile_bin(ii,ib,iparm))
3173 #endif
3174       endif
3175       if (nfile_asc(ii,ib,iparm).gt.0) then
3176         call card_concat(controlcard,.false.)
3177         call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
3178          maxfile_prot,nfile_asc(ii,ib,iparm))
3179 #ifdef DEBUG
3180         write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
3181         write(iout,*) (protfiles(i,2,ii,ib,iparm),&
3182           i=1,nfile_asc(ii,ib,iparm))
3183 #endif
3184       else if (nfile_cx(ii,ib,iparm).gt.0) then
3185         call card_concat(controlcard,.false.)
3186         call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
3187          maxfile_prot,nfile_cx(ii,ib,iparm))
3188 #ifdef DEBUG
3189         write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
3190         write(iout,*) (protfiles(i,2,ii,ib,iparm),&
3191          i=1,nfile_cx(ii,ib,iparm))
3192 #endif
3193       endif
3194       call flush(iout)
3195       enddo
3196       enddo
3197
3198       enddo
3199       return
3200       end subroutine read_protein_data
3201 !-------------------------------------------------------------------------------
3202       subroutine readsss(rekord,lancuch,wartosc,default)
3203 !      implicit none
3204       character*(*) :: rekord,lancuch,wartosc,default
3205       character(len=80) :: aux
3206       integer :: lenlan,lenrec,iread,ireade
3207 !el      external ilen
3208 !el      logical iblnk
3209 !el      external iblnk
3210       lenlan=ilen(lancuch)
3211       lenrec=ilen(rekord)
3212       iread=index(rekord,lancuch(:lenlan)//"=")
3213 !      print *,"rekord",rekord," lancuch",lancuch
3214 !      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
3215       if (iread.eq.0) then
3216         wartosc=default
3217         return
3218       endif
3219       iread=iread+lenlan+1
3220 !      print *,"iread",iread
3221 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3222       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
3223         iread=iread+1
3224 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3225       enddo
3226 !      print *,"iread",iread
3227       if (iread.gt.lenrec) then
3228          wartosc=default
3229         return
3230       endif
3231       ireade=iread+1
3232 !      print *,"ireade",ireade
3233       do while (ireade.lt.lenrec .and. &
3234          .not.iblnk(rekord(ireade:ireade)))
3235         ireade=ireade+1
3236       enddo
3237       wartosc=rekord(iread:ireade)
3238       return
3239       end subroutine readsss
3240 !----------------------------------------------------------------------------
3241       subroutine multreads(rekord,lancuch,tablica,dim,default)
3242 !      implicit none
3243       integer :: dim,i
3244       character*(*) rekord,lancuch,tablica(dim),default
3245       character(len=80) :: aux
3246       integer :: lenlan,lenrec,iread,ireade
3247 !el      external ilen
3248 !el      logical iblnk
3249 !el      external iblnk
3250       do i=1,dim
3251         tablica(i)=default
3252       enddo
3253       lenlan=ilen(lancuch)
3254       lenrec=ilen(rekord)
3255       iread=index(rekord,lancuch(:lenlan)//"=")
3256 !      print *,"rekord",rekord," lancuch",lancuch
3257 !      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
3258       if (iread.eq.0) return
3259       iread=iread+lenlan+1
3260       do i=1,dim
3261 !      print *,"iread",iread
3262 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3263       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
3264         iread=iread+1
3265 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3266       enddo
3267 !      print *,"iread",iread
3268       if (iread.gt.lenrec) return
3269       ireade=iread+1
3270 !      print *,"ireade",ireade
3271       do while (ireade.lt.lenrec .and. &
3272          .not.iblnk(rekord(ireade:ireade)))
3273         ireade=ireade+1
3274       enddo
3275       tablica(i)=rekord(iread:ireade)
3276       iread=ireade+1
3277       enddo
3278       end subroutine multreads
3279 !----------------------------------------------------------------------------
3280       subroutine split_string(rekord,tablica,dim,nsub)
3281 !      implicit none
3282       integer :: dim,nsub,i,ii,ll,kk
3283       character*(*) tablica(dim)
3284       character*(*) rekord
3285 !el      integer ilen
3286 !el      external ilen
3287       do i=1,dim
3288         tablica(i)=" "
3289       enddo
3290       ii=1
3291       ll = ilen(rekord)
3292       nsub=0
3293       do i=1,dim
3294 ! Find the start of term name
3295         kk = 0
3296         do while (ii.le.ll .and. rekord(ii:ii).eq." ")
3297           ii = ii+1
3298         enddo
3299 ! Parse the name into TABLICA(i) until blank found
3300         do while (ii.le.ll .and. rekord(ii:ii).ne." ")
3301           kk = kk+1
3302           tablica(i)(kk:kk)=rekord(ii:ii)
3303           ii = ii+1
3304         enddo
3305         if (kk.gt.0) nsub=nsub+1
3306         if (ii.gt.ll) return
3307       enddo
3308       return
3309       end subroutine split_string
3310 !--------------------------------------------------------------------------------
3311 ! readrtns_compar.F
3312 !--------------------------------------------------------------------------------
3313       subroutine read_compar
3314 !
3315 ! Read molecular data
3316 !
3317       use conform_compar, only:alloc_compar_arrays
3318       use control_data, only:pdbref
3319       use geometry_data, only:deg2rad,rad2deg
3320 !      implicit none
3321 !      include 'DIMENSIONS'
3322 !      include 'DIMENSIONS.ZSCOPT'
3323 !      include 'DIMENSIONS.COMPAR'
3324 !      include 'DIMENSIONS.FREE'
3325 !      include 'COMMON.IOUNITS'
3326 !      include 'COMMON.TIME1'
3327 !      include 'COMMON.SBRIDGE'
3328 !      include 'COMMON.CONTROL'
3329 !      include 'COMMON.COMPAR'
3330 !      include 'COMMON.CHAIN'
3331 !      include 'COMMON.HEADER'
3332 !      include 'COMMON.GEO'
3333 !      include 'COMMON.FREE'
3334       character(len=320) :: controlcard !,ucase
3335       character(len=64) :: wfile
3336 !el      integer ilen
3337 !el      external ilen
3338       integer :: i,j,k
3339 !elwrite(iout,*)"jestesmy w read_compar"
3340       call card_concat(controlcard,.true.)
3341       pdbref=(index(controlcard,'PDBREF').gt.0)
3342       call reada(controlcard,'CUTOFF_UP',rmscut_base_up,4.0d0)
3343       call reada(controlcard,'CUTOFF_LOW',rmscut_base_low,3.0d0)
3344       call reada(controlcard,'RMSUP_LIM',rmsup_lim,4.0d0)
3345       call reada(controlcard,'RMSUPUP_LIM',rmsupup_lim,7.5d0)
3346       verbose = index(controlcard,"VERBOSE").gt.0
3347       lgrp=index(controlcard,"STATIN").gt.0
3348       lgrp_out=index(controlcard,"STATOUT").gt.0
3349       merge_helices=index(controlcard,"DONT_MERGE_HELICES").eq.0
3350       binary = index(controlcard,"BINARY").gt.0
3351       rmscut_base_up=rmscut_base_up/50
3352       rmscut_base_low=rmscut_base_low/50
3353       call reada(controlcard,"FRAC_SEC",frac_sec,0.66666666d0)
3354       call readi(controlcard,'NLEVEL',nlevel,1)
3355       if (nlevel.lt.0) then
3356         allocate(nfrag(2))
3357         call alloc_compar_arrays(maxfrag,1)
3358         goto 121
3359       else
3360         allocate(nfrag(nlevel))
3361       endif
3362 ! Read the data pertaining to elementary fragments (level 1)
3363       call readi(controlcard,'NFRAG',nfrag(1),0)
3364       write(iout,*)"nfrag(1)",nfrag(1)
3365       call alloc_compar_arrays(nfrag(1),nlevel)
3366       do j=1,nfrag(1)
3367         call card_concat(controlcard,.true.)
3368         write (iout,*) controlcard(:ilen(controlcard))
3369         call readi(controlcard,'NPIECE',npiece(j,1),0)
3370         call readi(controlcard,'N_SHIFT1',n_shift(1,j,1),0)
3371         call readi(controlcard,'N_SHIFT2',n_shift(2,j,1),0)
3372         call reada(controlcard,'ANGCUT',ang_cut(j),50.0d0)
3373         call reada(controlcard,'MAXANG',ang_cut1(j),360.0d0)
3374         call reada(controlcard,'FRAC_MIN',frac_min(j),0.666666d0)
3375         call reada(controlcard,'NC_FRAC',nc_fragm(j,1),0.5d0)
3376         call readi(controlcard,'NC_REQ',nc_req_setf(j,1),0)
3377         call readi(controlcard,'RMS',irms(j,1),0)
3378         call readi(controlcard,'LOCAL',iloc(j),1)
3379         call readi(controlcard,'ELCONT',ielecont(j,1),1)
3380         if (ielecont(j,1).eq.0) then
3381           call readi(controlcard,'SCCONT',isccont(j,1),1)
3382         endif
3383         ang_cut(j)=ang_cut(j)*deg2rad
3384         ang_cut1(j)=ang_cut1(j)*deg2rad
3385         do k=1,npiece(j,1)
3386           call card_concat(controlcard,.true.)
3387           call readi(controlcard,'IFRAG1',ifrag(1,k,j),0)
3388           call readi(controlcard,'IFRAG2',ifrag(2,k,j),0)
3389         enddo
3390         write(iout,*)"j",j," npiece",npiece(j,1)," ifrag",&
3391           (ifrag(1,k,j),ifrag(2,k,j),&
3392          k=1,npiece(j,1))," ang_cut",ang_cut(j)*rad2deg,&
3393           " ang_cut1",ang_cut1(j)*rad2deg
3394         write(iout,*)"n_shift",n_shift(1,j,1),n_shift(2,j,1)
3395         write(iout,*)"nc_frac",nc_fragm(j,1)," nc_req",nc_req_setf(j,1)
3396         write(iout,*)"irms",irms(j,1)," ielecont",ielecont(j,1),&
3397           " ilocal",iloc(j)," isccont",isccont(j,1)
3398       enddo
3399 ! Read data pertaning to higher levels
3400       do i=2,nlevel
3401         call card_concat(controlcard,.true.)
3402         call readi(controlcard,'NFRAG',NFRAG(i),0)
3403         write (iout,*) "i",i," nfrag",nfrag(i)
3404         do j=1,nfrag(i)
3405           call card_concat(controlcard,.true.)
3406           if (i.eq.2) then
3407             call readi(controlcard,'ELCONT',ielecont(j,i),0)
3408             if (ielecont(j,i).eq.0) then
3409               call readi(controlcard,'SCCONT',isccont(j,i),1)
3410             endif
3411             call readi(controlcard,'RMS',irms(j,i),0)
3412           else
3413             ielecont(j,i)=0
3414             isccont(j,i)=0
3415             irms(j,i)=1
3416           endif
3417           call readi(controlcard,'NPIECE',npiece(j,i),0)
3418           call readi(controlcard,'N_SHIFT1',n_shift(1,j,i),0)
3419           call readi(controlcard,'N_SHIFT2',n_shift(2,j,i),0)
3420           call multreadi(controlcard,'IPIECE',ipiece(1,j,i),&
3421             npiece(j,i),0)
3422           call reada(controlcard,'NC_FRAC',nc_fragm(j,i),0.5d0)
3423           call readi(controlcard,'NC_REQ',nc_req_setf(j,i),0)
3424           write(iout,*) "j",j," npiece",npiece(j,i)," n_shift",&
3425             n_shift(1,j,i),n_shift(2,j,i)," ielecont",ielecont(j,i),&
3426             " isccont",isccont(j,i)," irms",irms(j,i)
3427           write(iout,*) "ipiece",(ipiece(k,j,i),k=1,npiece(j,i))
3428           write(iout,*)"n_shift",n_shift(1,j,i),n_shift(2,j,i)
3429           write(iout,*)"nc_frac",nc_fragm(j,i),&
3430            " nc_req",nc_req_setf(j,i)
3431         enddo
3432       enddo
3433       if (binary) write (iout,*) "Classes written in binary format."
3434       return
3435   121 continue
3436       call reada(controlcard,'ANGCUT_HEL',angcut_hel,50.0d0)
3437       call reada(controlcard,'MAXANG_HEL',angcut1_hel,60.0d0)
3438       call reada(controlcard,'ANGCUT_BET',angcut_bet,90.0d0)
3439       call reada(controlcard,'MAXANG_BET',angcut1_bet,360.0d0)
3440       call reada(controlcard,'ANGCUT_STRAND',angcut_strand,90.0d0)
3441       call reada(controlcard,'MAXANG_STRAND',angcut1_strand,60.0d0)
3442       call reada(controlcard,'FRAC_MIN',frac_min_set,0.666666d0)
3443       call reada(controlcard,'NC_FRAC_HEL',ncfrac_hel,0.5d0)
3444       call readi(controlcard,'NC_REQ_HEL',ncreq_hel,0)
3445       call reada(controlcard,'NC_FRAC_BET',ncfrac_bet,0.5d0)
3446       call reada(controlcard,'NC_FRAC_PAIR',ncfrac_pair,0.3d0)
3447       call readi(controlcard,'NC_REQ_BET',ncreq_bet,0)
3448       call readi(controlcard,'NC_REQ_PAIR',ncreq_pair,0)
3449       call readi(controlcard,'NSHIFT_HEL',nshift_hel,3)
3450       call readi(controlcard,'NSHIFT_BET',nshift_bet,3)
3451       call readi(controlcard,'NSHIFT_STRAND',nshift_strand,3)
3452       call readi(controlcard,'NSHIFT_PAIR',nshift_pair,3)
3453       call readi(controlcard,'RMS_SINGLE',irms_single,0)
3454       call readi(controlcard,'CONT_SINGLE',icont_single,1)
3455       call readi(controlcard,'LOCAL_SINGLE',iloc_single,1)
3456       call readi(controlcard,'RMS_PAIR',irms_pair,0)
3457       call readi(controlcard,'CONT_PAIR',icont_pair,1)
3458       call readi(controlcard,'SPLIT_BET',isplit_bet,0)
3459       angcut_hel=angcut_hel*deg2rad
3460       angcut1_hel=angcut1_hel*deg2rad
3461       angcut_bet=angcut_bet*deg2rad
3462       angcut1_bet=angcut1_bet*deg2rad
3463       angcut_strand=angcut_strand*deg2rad
3464       angcut1_strand=angcut1_strand*deg2rad
3465       write (iout,*) "Automatic detection of structural elements"
3466       write (iout,*) 'NC_FRAC_HEL',ncfrac_hel,' NC_REQ_HEL',ncreq_hel,&
3467                      ' NC_FRAC_BET',ncfrac_bet,' NC_REQ_BET',ncreq_bet,&
3468                  ' RMS_SINGLE',irms_single,' CONT_SINGLE',icont_single,&
3469                  ' NC_FRAC_PAIR',ncfrac_pair,' NC_REQ_PAIR',ncreq_pair,&
3470         ' RMS_PAIR',irms_pair,' CONT_PAIR',icont_pair,&
3471         ' SPLIT_BET',isplit_bet
3472       write (iout,*) 'NSHIFT_HEL',nshift_hel,' NSHIFT_BET',nshift_bet,&
3473         ' NSHIFT_STRAND',nshift_strand,' NSHIFT_PAIR',nshift_pair
3474       write (iout,*) 'ANGCUT_HEL',angcut_hel*rad2deg,&
3475         ' MAXANG_HEL',angcut1_hel*rad2deg
3476       write (iout,*) 'ANGCUT_BET',angcut_bet*rad2deg,&
3477                      ' MAXANG_BET',angcut1_bet*rad2deg
3478       write (iout,*) 'ANGCUT_STRAND',angcut_strand*rad2deg,&
3479                      ' MAXANG_STRAND',angcut1_strand*rad2deg
3480       write (iout,*) 'FRAC_MIN',frac_min_set
3481       return
3482       end subroutine read_compar
3483 !--------------------------------------------------------------------------------
3484 ! read_ref_str.F
3485 !--------------------------------------------------------------------------------
3486       subroutine read_ref_structure(*)
3487 !
3488 ! Read the reference structure from the PDB file or from a PDB file or in the form of the dihedral
3489 ! angles.
3490 !
3491       use control_data, only:pdbref
3492       use geometry_data, only:nres,cref,c,dc,nsup,dc_norm,nend_sup,&
3493                               nstart_sup,nstart_seq,nperm,nres0
3494       use energy_data, only:nct,nnt,icont_ref,ncont_ref,itype,molnum
3495       use compare, only:seq_comp !,contact,elecont
3496       use geometry, only:chainbuild,dist
3497       use io_config, only:readpdb
3498 !
3499       use conform_compar, only:contact,elecont
3500 !      implicit none
3501 !      include 'DIMENSIONS'
3502 !      include 'DIMENSIONS.ZSCOPT'
3503 !      include 'DIMENSIONS.COMPAR'
3504 !      include 'COMMON.IOUNITS'
3505 !      include 'COMMON.GEO'
3506 !      include 'COMMON.VAR'
3507 !      include 'COMMON.INTERACT'
3508 !      include 'COMMON.LOCAL'
3509 !      include 'COMMON.NAMES'
3510 !      include 'COMMON.CHAIN'
3511 !      include 'COMMON.FFIELD'
3512 !      include 'COMMON.SBRIDGE'
3513 !      include 'COMMON.HEADER'
3514 !      include 'COMMON.CONTROL'
3515 !      include 'COMMON.CONTACTS1'
3516 !      include 'COMMON.PEPTCONT'
3517 !      include 'COMMON.TIME1'
3518 !      include 'COMMON.COMPAR'
3519       character(len=4) :: sequence(nres)
3520 !el      integer rescode
3521 !el      real(kind=8) :: x(maxvar)
3522       integer :: itype_pdb(nres,5)
3523 !el      logical seq_comp
3524       integer :: i,j,k,nres_pdb,iaux,mnum
3525       real(kind=8) :: ddsc !el,dist
3526       integer :: kkk !,ilen
3527 !el      external ilen
3528 !
3529       nres0=nres
3530       write (iout,*) "pdbref",pdbref
3531       if (pdbref) then
3532         read(inp,'(a)') pdbfile
3533         write (iout,'(2a,1h.)') 'PDB data will be read from file ',&
3534           pdbfile(:ilen(pdbfile))
3535         open(ipdbin,file=pdbfile,status='old',err=33)
3536         goto 34
3537   33    write (iout,'(a)') 'Error opening PDB file.'
3538         return 1
3539   34    continue
3540         do i=1,nres
3541           mnum=molnum(i)
3542           itype_pdb(i,mnum)=itype(i,mnum)
3543         enddo
3544
3545         call readpdb
3546
3547         do i=1,nres
3548           iaux=itype_pdb(i,mnum)
3549           itype_pdb(i,mnum)=itype(i,mnum)
3550           itype(i,mnum)=iaux
3551         enddo
3552         close (ipdbin)
3553         do kkk=1,nperm
3554         nres_pdb=nres
3555         nres=nres0
3556         nstart_seq=nnt
3557         if (nsup.le.(nct-nnt+1)) then
3558           do i=0,nct-nnt+1-nsup
3559             if (seq_comp(itype(nnt+i,molnum(nnt+i)),itype_pdb(nstart_sup,molnum(nstart_sup)),&
3560               nsup)) then
3561               do j=nnt+nsup-1,nnt,-1
3562                 do k=1,3
3563                   cref(k,nres+j+i,kkk)=cref(k,nres_pdb+j,kkk)
3564                 enddo
3565               enddo
3566               do j=nnt+nsup-1,nnt,-1
3567                 do k=1,3
3568                   cref(k,j+i,kkk)=cref(k,j,kkk)
3569                 enddo
3570                 write(iout,*) "J",j,"J+I",j+i
3571                 phi_ref(j+i)=phi_ref(j)
3572                 theta_ref(j+i)=theta_ref(j)
3573                 alph_ref(j+i)=alph_ref(j)
3574                 omeg_ref(j+i)=omeg_ref(j)
3575               enddo
3576 #ifdef DEBUG
3577               do j=nnt,nct
3578                 write (iout,'(i5,3f10.5,5x,3f10.5)') &
3579                   j,(cref(k,j,kkk),k=1,3),(cref(k,j+nres,kkk),k=1,3)
3580               enddo
3581 #endif
3582               nstart_seq=nnt+i
3583               nstart_sup=nnt+i
3584               goto 111
3585             endif
3586           enddo
3587           write (iout,'(a)') &
3588                   'Error - sequences to be superposed do not match.'
3589           return 1
3590         else
3591           do i=0,nsup-(nct-nnt+1)
3592             if (seq_comp(itype(nnt,molnum(nnt)),itype_pdb(nstart_sup+i,molnum(nstart_sup+i)),&
3593               nct-nnt+1)) &
3594             then
3595               nstart_sup=nstart_sup+i
3596               nsup=nct-nnt+1
3597               goto 111
3598             endif
3599           enddo
3600           write (iout,'(a)') &
3601                   'Error - sequences to be superposed do not match.'
3602         endif
3603         enddo
3604   111   continue
3605         write (iout,'(a,i5)') &
3606          'Experimental structure begins at residue',nstart_seq
3607       else
3608         call read_angles(inp,*38)
3609         goto 39
3610    38   write (iout,'(a)') 'Error reading reference structure.'
3611         return 1
3612    39   call chainbuild
3613         kkk=1
3614         nstart_sup=nnt
3615         nstart_seq=nnt
3616         nsup=nct-nnt+1
3617         do i=1,2*nres
3618           do j=1,3
3619             cref(j,i,kkk)=c(j,i)
3620           enddo
3621         enddo
3622       endif
3623       nend_sup=nstart_sup+nsup-1
3624       do i=1,2*nres
3625         do j=1,3
3626           c(j,i)=cref(j,i,kkk)
3627         enddo
3628       enddo
3629       do i=1,nres
3630         mnum=molnum(i)
3631         do j=1,3
3632           dc(j,nres+i)=cref(j,nres+i,kkk)-cref(j,i,kkk)
3633         enddo
3634         if (itype(i,mnum).ne.10) then
3635           ddsc = dist(i,nres+i)
3636           do j=1,3
3637             dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
3638           enddo
3639         else
3640           do j=1,3
3641             dc_norm(j,nres+i)=0.0d0
3642           enddo
3643         endif
3644 !        write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
3645 !         " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
3646 !         dc_norm(3,nres+i)**2
3647         do j=1,3
3648           dc(j,i)=c(j,i+1)-c(j,i)
3649         enddo
3650         ddsc = dist(i,i+1)
3651         do j=1,3
3652           dc_norm(j,i)=dc(j,i)/ddsc
3653         enddo
3654       enddo
3655 !      print *,"Calling contact"
3656       call contact(.true.,ncont_ref,icont_ref(1,1),&
3657         nstart_sup,nend_sup)
3658 !      print *,"Calling elecont"
3659       call elecont(.true.,ncont_pept_ref,&
3660          icont_pept_ref(1,1),&
3661          nstart_sup,nend_sup)
3662        write (iout,'(a,i3,a,i3,a,i3,a)') &
3663           'Number of residues to be superposed:',nsup,&
3664           ' (from residue',nstart_sup,' to residue',&
3665           nend_sup,').'
3666       return
3667       end subroutine read_ref_structure
3668 !--------------------------------------------------------------------------------
3669 ! geomout.F
3670 !--------------------------------------------------------------------------------
3671       subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev)
3672
3673       use geometry_data, only:nres,c
3674       use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype,molnum
3675 !      implicit real*8 (a-h,o-z)
3676 !      include 'DIMENSIONS'
3677 !      include 'DIMENSIONS.ZSCOPT'
3678 !      include 'COMMON.CHAIN'
3679 !      include 'COMMON.INTERACT'
3680 !      include 'COMMON.NAMES'
3681 !      include 'COMMON.IOUNITS'
3682 !      include 'COMMON.HEADER'
3683 !      include 'COMMON.SBRIDGE'
3684       character(len=50) :: tytul
3685       character(len=1),dimension(10) :: chainid=reshape((/'A','B','C',&
3686                       'D','E','F','G','H','I','J'/),shape(chainid))
3687       integer,dimension(nres) :: ica !(maxres)
3688       real(kind=8) :: temp,efree,etot,entropy,rmsdev
3689       integer :: ii,i,j,iti,ires,iatom,ichain,mnum
3690       write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')&
3691         ii,temp,rmsdev
3692       write (ipdb,'("REMARK DIMENSIONLESS FREE ENERGY",1pe15.5)') &
3693         efree
3694       write (ipdb,'("REMARK ENERGY",1pe15.5," ENTROPY",1pe15.5)') &
3695         etot,entropy
3696       iatom=0
3697       ichain=1
3698       ires=0
3699       do i=nnt,nct
3700         mnum=molnum(i)
3701         iti=itype(i,mnum)
3702         if (iti.eq.ntyp1) then
3703           ichain=ichain+1
3704           ires=0
3705           write (ipdb,'(a)') 'TER'
3706         else
3707         ires=ires+1
3708         iatom=iatom+1
3709         ica(i)=iatom
3710         write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
3711            ires,(c(j,i),j=1,3)
3712         if (iti.ne.10) then
3713           iatom=iatom+1
3714           write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
3715             ires,(c(j,nres+i),j=1,3)
3716         endif
3717         endif
3718       enddo
3719       write (ipdb,'(a)') 'TER'
3720       do i=nnt,nct-1
3721         mnum=molnum(i)
3722         if (itype(i,mnum).eq.ntyp1) cycle
3723         if (itype(i,mnum).eq.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
3724           write (ipdb,30) ica(i),ica(i+1)
3725         else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
3726           write (ipdb,30) ica(i),ica(i+1),ica(i)+1
3727         else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).eq.ntyp1_molec(mnum)) then
3728           write (ipdb,30) ica(i),ica(i)+1
3729         endif
3730       enddo
3731       if (itype(nct,molnum(nct)).ne.10) then
3732         write (ipdb,30) ica(nct),ica(nct)+1
3733       endif
3734       do i=1,nss
3735         write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
3736       enddo
3737       write (ipdb,'(a6)') 'ENDMDL'
3738   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
3739   20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
3740   30  FORMAT ('CONECT',8I5)
3741       return
3742       end subroutine pdboutW
3743 #endif
3744 !------------------------------------------------------------------------------
3745       end module io_wham
3746 !-----------------------------------------------------------------------------
3747 !-----------------------------------------------------------------------------
3748