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       call mygetenv('IONPAR_NUCL',ionnuclname)
 102       open (iionnucl,file=ionnuclname,status='old')
 103       call mygetenv('IONPAR_TRAN',iontranname)
 104       open (iiontran,file=iontranname,status='old',action='read')
 105
 106
 107 #ifndef OLDSCP
 108 !
 109 ! 8/9/01 In the newest version SCp interaction constants are read from a file
 110 ! Use -DOLDSCP to use hard-coded constants instead.
 111 !
 112       call mygetenv('SCPPAR',scpname)
 113       open (iscpp,file=scpname,status='old')
 114 #endif
 115 #ifdef MPL
 116       if (MyID.eq.BossID) then
 117       MyRank = MyID/fgProcs
 118 #endif
 119 #ifdef MPI
 120       print *,'OpenUnits: processor',MyRank
 121       call numstr(MyRank,liczba)
 122       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//liczba
 123 #else
 124       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
 125 #endif
 126       open(iout,file=outname,status='unknown')
 127       write (iout,'(80(1h-))')
 128       write (iout,'(30x,a)') "FILE ASSIGNMENT"
 129       write (iout,'(80(1h-))')
 130       write (iout,*) "Input file                      : ",&
 131         prefix(:ilen(prefix))//'.inp'
 132       write (iout,*) "Output file                     : ",&
 133         outname(:ilen(outname))
 134       write (iout,*)
 135       write (iout,*) "Sidechain potential file        : ",&
 136         sidename(:ilen(sidename))
 137 #ifndef OLDSCP
 138       write (iout,*) "SCp potential file              : ",&
 139         scpname(:ilen(scpname))
 140 #endif
 141       write (iout,*) "Electrostatic potential file    : ",&
 142         elename(:ilen(elename))
 143       write (iout,*) "Cumulant coefficient file       : ",&
 144         fouriername(:ilen(fouriername))
 145       write (iout,*) "Torsional parameter file        : ",&
 146         torname(:ilen(torname))
 147       write (iout,*) "Double torsional parameter file : ",&
 148         tordname(:ilen(tordname))
 149       write (iout,*) "Backbone-rotamer parameter file : ",&
 150         sccorname(:ilen(sccorname))
 151       write (iout,*) "Bond & inertia constant file    : ",&
 152         bondname(:ilen(bondname))
 153       write (iout,*) "Bending parameter file          : ",&
 154         thetname(:ilen(thetname))
 155       write (iout,*) "Rotamer parameter file          : ",&
 156         rotname(:ilen(rotname))
 157       write (iout,'(80(1h-))')
 158       write (iout,*)
 159       return
 160       end subroutine openunits
 161 !-----------------------------------------------------------------------------
 162 ! molread_zs.F
 163 !-----------------------------------------------------------------------------
 164       subroutine molread(*)
 165 !
 166 ! Read molecular data.
 167 !
 168       use energy_data
 169       use geometry_data, only:nres,deg2rad,c,dc,nres_molec,crefjlee,cref
 170       use control_data, only:iscode,dyn_ss,pdbref,indpdb
 171       use io_base, only:rescode
 172       use control, only:setup_var,init_int_table,hpb_partition
 173       use geometry, only:alloc_geo_arrays
 174       use energy, only:alloc_ener_arrays
 175 !      implicit real*8 (a-h,o-z)
 176 !      include 'DIMENSIONS'
 177 !      include 'DIMENSIONS.ZSCOPT'
 178 !      include 'COMMON.IOUNITS'
 179 !      include 'COMMON.GEO'
 180 !      include 'COMMON.VAR'
 181 !      include 'COMMON.INTERACT'
 182 !      include 'COMMON.LOCAL'
 183 !      include 'COMMON.NAMES'
 184 !      include 'COMMON.CHAIN'
 185 !      include 'COMMON.FFIELD'
 186 !      include 'COMMON.SBRIDGE'
 187 !      include 'COMMON.TORCNSTR'
 188 !      include 'COMMON.CONTROL'
 189       character(len=4),dimension(:),allocatable :: sequence !(nres)
 190 !el      integer :: rescode
 191 !el      real(kind=8) :: x(maxvar)
 192       character(len=320) :: controlcard !,ucase
 193       integer,dimension(nres,5) :: itype_pdb !(maxres)
 194       integer :: i,j,i1,i2,it1,it2,mnum
 195       real(kind=8) :: scalscp
 196 !el      logical :: seq_comp
 197       write(iout,*) "START MOLREAD"
 198       call flush(iout)
 199       do i=1,5
 200        nres_molec(i)=0
 201        print *,"nres_molec, initial",nres_molec(i)
 202       enddo
 203
 204       call card_concat(controlcard,.true.)
 205       call reada(controlcard,'SCAL14',scal14,0.4d0)
 206       call reada(controlcard,'SCALSCP',scalscp,1.0d0)
 207       call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
 208       call reada(controlcard,'TEMP0',temp0,300.0d0) !el
 209       call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
 210       r0_corr=cutoff_corr-delt_corr
 211       call readi(controlcard,"NRES",nres_molec(1),0)
 212       call readi(controlcard,"NRES_NUCL",nres_molec(2),0)
 213       call readi(controlcard,"NRES_CAT",nres_molec(5),0)
 214       nres=0
 215       do i=1,5
 216        nres=nres_molec(i)+nres
 217       enddo
 218 !      allocate(sequence(nres+1))
 219 !el znamy juz ilosc reszt wiec mozna zaalokowac tablice do liczenia enerii
 220       call alloc_geo_arrays
 221       call alloc_ener_arrays
 222 ! alokacja dodatkowych tablic, ktore w unresie byly alokowanie w locie
 223 !----------------------------
 224       allocate(c(3,2*nres+2))
 225       allocate(dc(3,0:2*nres+2))
 226       allocate(itype(nres+2,5))
 227       allocate(itel(nres+2))
 228       if (.not. allocated(molnum)) allocate(molnum(nres+2))
 229 !
 230 ! Zero out tableis.
 231       do i=1,2*nres+2
 232         do j=1,3
 233           c(j,i)=0.0D0
 234           dc(j,i)=0.0D0
 235         enddo
 236       enddo
 237       do i=1,nres+2
 238         molnum(i)=0
 239         do j=1,5
 240         itype(i,j)=0
 241         enddo
 242         itel(i)=0
 243       enddo
 244 !--------------------------
 245 !
 246       allocate(sequence(nres+1))
 247       iscode=index(controlcard,"ONE_LETTER")
 248       if (nres.le.0) then
 249         write (iout,*) "Error: no residues in molecule"
 250         return 1
 251       endif
 252       if (nres.gt.maxres) then
 253         write (iout,*) "Error: too many residues",nres,maxres
 254       endif
 255       write(iout,*) 'nres=',nres
 256 ! HERE F**** CHANGE
 257 ! Read sequence of the protein
 258       if (iscode.gt.0) then
 259         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 260       else
 261         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 262       endif
 263 ! Convert sequence to numeric code
 264       do i=1,nres_molec(1)
 265         mnum=1
 266         molnum(i)=1
 267         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
 268       enddo
 269       do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
 270         mnum=2
 271         molnum(i)=2
 272         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
 273       enddo
 274       do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
 275         mnum=5
 276         molnum(i)=5
 277         itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
 278       enddo
 279
 280       write (iout,*) "Numeric code:"
 281       write (iout,'(20i4)') (itype(i,molnum(i)),i=1,nres)
 282       do i=1,nres-1
 283         mnum=molnum(i)
 284 #ifdef PROCOR
 285         if (itype(i,mnum).eq.ntyp1_molec(mnum) .or. itype(i+1,mnum).eq.ntyp1_molec(mnum)) then
 286 #else
 287         if (itype(i,mnum).eq.ntyp1_molec(mnum)) then
 288 #endif
 289           itel(i)=0
 290 #ifdef PROCOR
 291         else if (iabs(itype(i+1,mnum)).ne.20) then
 292 #else
 293         else if (iabs(itype(i,mnum)).ne.20) then
 294 #endif
 295           itel(i)=1
 296         else
 297           itel(i)=2
 298         endif
 299       enddo
 300        write (iout,*) "ITEL"
 301        do i=1,nres-1
 302          mnum=molnum(i)
 303          write (iout,*) i,itype(i,mnum),itel(i)
 304        enddo
 305       write(iout,*)
 306       call read_bridge
 307
 308       if (with_dihed_constr) then
 309
 310       read (inp,*) ndih_constr
 311       if (ndih_constr.gt.0) then
 312         read (inp,*) ftors
 313         write (iout,*) 'FTORS',ftors
 314         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 315         write (iout,*) &
 316          'There are',ndih_constr,' constraints on phi angles.'
 317         do i=1,ndih_constr
 318           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 319         enddo
 320         do i=1,ndih_constr
 321           phi0(i)=deg2rad*phi0(i)
 322           drange(i)=deg2rad*drange(i)
 323         enddo
 324       endif
 325
 326       endif
 327       nnt=1
 328       nct=nres
 329       if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
 330       if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
 331       write(iout,*) 'NNT=',NNT,' NCT=',NCT
 332       if (constr_homology.gt.0) then
 333 !c       write (iout,*) "About to call read_constr_homology"
 334 !c       call flush(iout)
 335         call read_constr_homology
 336 !c       write (iout,*) "Exit read_constr_homology"
 337 !c       call flush(iout)
 338         if (indpdb.gt.0 .or. pdbref) then
 339           do i=1,2*nres
 340             do j=1,3
 341               c(j,i)=crefjlee(j,i)
 342               cref(j,i,1)=crefjlee(j,i)
 343             enddo
 344           enddo
 345         endif
 346         endif
 347 #ifdef DEBUG
 348         write (iout,*) "Array C"
 349         do i=1,nres
 350           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
 351      &      (c(j,i+nres),j=1,3)
 352         enddo
 353         write (iout,*) "Array Cref"
 354         do i=1,nres
 355           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),
 356      &      (cref(j,i+nres,1),j=1,3)
 357         enddo
 358 #endif
 359       call setup_var
 360       call init_int_table
 361       if (ns.gt.0) then
 362         write (iout,'(/a,i3,a)') 'The chain contains',ns,&
 363         ' disulfide-bridging cysteines.'
 364         write (iout,'(20i4)') (iss(i),i=1,ns)
 365         write (iout,'(/a/)') 'Pre-formed links are:'
 366         do i=1,nss
 367           mnum=molnum(i)
 368           i1=ihpb(i)-nres
 369           i2=jhpb(i)-nres
 370           it1=itype(i1,mnum)
 371           it2=itype(i2,mnum)
 372          write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
 373           restyp(it1,molnum(i1)),'(',i1,') -- ',restyp(it2,molnum(i2)),'(',i2,')',&
 374           dhpb(i),ebr,forcon(i)
 375         enddo
 376       endif
 377       if (ns.gt.0.and.dyn_ss) then
 378           do i=nss+1,nhpb
 379             ihpb(i-nss)=ihpb(i)
 380             jhpb(i-nss)=jhpb(i)
 381             forcon(i-nss)=forcon(i)
 382             dhpb(i-nss)=dhpb(i)
 383           enddo
 384           nhpb=nhpb-nss
 385           nss=0
 386           call hpb_partition
 387           do i=1,ns
 388             dyn_ss_mask(iss(i))=.true.
 389           enddo
 390       endif
 391       write (iout,'(a)')
 392       return
 393       end subroutine molread
 394 !-----------------------------------------------------------------------------
 395 ! parmread.F
 396 !-----------------------------------------------------------------------------
 397       subroutine parmread(iparm,*)
 398 #else
 399       subroutine parmread
 400 #endif
 401 !
 402 ! Read the parameters of the probability distributions of the virtual-bond
 403 ! valence angles and the side chains and energy parameters.
 404 !
 405       use wham_data
 406
 407       use geometry_data
 408       use energy_data
 409       use control_data, only: maxterm,maxlor,maxterm_sccor,& !maxtor
 410           maxtermd_1,maxtermd_2,tor_mode,scelemode !,maxthetyp,maxthetyp1
 411       use MD_data
 412 !el      use MPI_data
 413 !el      use map_data
 414       use io_config, only: printmat
 415       use control, only: getenv_loc
 416
 417 #ifdef MPI
 418       use MPI_data
 419       include "mpif.h"
 420       integer :: IERROR
 421 #endif
 422 !      implicit real*8 (a-h,o-z)
 423 !      include 'DIMENSIONS'
 424 !      include 'DIMENSIONS.ZSCOPT'
 425 !      include 'DIMENSIONS.FREE'
 426 !      include 'COMMON.IOUNITS'
 427 !      include 'COMMON.CHAIN'
 428 !      include 'COMMON.INTERACT'
 429 !      include 'COMMON.GEO'
 430 !      include 'COMMON.LOCAL'
 431 !      include 'COMMON.TORSION'
 432 !      include 'COMMON.FFIELD'
 433 !      include 'COMMON.NAMES'
 434 !      include 'COMMON.SBRIDGE'
 435 !      include 'COMMON.WEIGHTS'
 436 !      include 'COMMON.ENEPS'
 437 !      include 'COMMON.SCCOR'
 438 !      include 'COMMON.SCROT'
 439 !      include 'COMMON.FREE'
 440       character(len=1) :: t1,t2,t3
 441       character(len=1) :: onelett(4) = (/"G","A","P","D"/)
 442       character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
 443       logical :: lprint,SPLIT_FOURIERTOR
 444       real(kind=8),dimension(3,3,maxlob) :: blower      !(3,3,maxlob)
 445       character(len=800) :: controlcard
 446       character(len=256) :: bondname_t,thetname_t,rotname_t,torname_t,&
 447         tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,&
 448         sccorname_t
 449 !el      integer ilen
 450 !el   external ilen
 451       character(len=16) :: key
 452       integer :: iparm,nkcctyp
 453 !el      real(kind=8) :: ip,mp
 454       real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,epsijlip,&
 455                 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,    &
 456                 epspeptube,epssctube,sigmapeptube,      &
 457                 sigmasctube,ssscale
 458
 459       real(kind=8) :: v0ij,v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,rjunk,&
 460                 res1
 461       integer :: i,j,ichir1,ichir2,k,l,m,kk,ii,mm,junk,lll,ll,llll,n,jj
 462       integer :: nlobi,iblock,maxinter,iscprol,ncatprotparm
 463       character*3 string
 464
 465 !
 466 ! Body
 467 !
 468 ! Set LPRINT=.TRUE. for debugging
 469       dwa16=2.0d0**(1.0d0/6.0d0)
 470       lprint=.false.
 471       itypro=20
 472 ! Assign virtual-bond length
 473       vbl=3.8D0
 474       vblinv=1.0D0/vbl
 475       vblinv2=vblinv*vblinv
 476       itime_mat=0
 477 #ifndef CLUSTER
 478       call card_concat(controlcard,.true.)
 479       wname(4)="WCORRH"
 480 !el
 481       call reada(controlcard,"D0CM",d0cm,3.78d0)
 482       call reada(controlcard,"AKCM",akcm,15.1d0)
 483       call reada(controlcard,"AKTH",akth,11.0d0)
 484       call reada(controlcard,"AKCT",akct,12.0d0)
 485       call reada(controlcard,"V1SS",v1ss,-1.08d0)
 486       call reada(controlcard,"V2SS",v2ss,7.61d0)
 487       call reada(controlcard,"V3SS",v3ss,13.7d0)
 488       call reada(controlcard,"EBR",ebr,-5.50D0)
 489       call reada(controlcard,"ATRISS",atriss,0.301D0)
 490       call reada(controlcard,"BTRISS",btriss,0.021D0)
 491       call reada(controlcard,"CTRISS",ctriss,1.001D0)
 492       call reada(controlcard,"DTRISS",dtriss,1.001D0)
 493       call reada(controlcard,"SSSCALE",ssscale,1.0D0)
 494       dyn_ss=(index(controlcard,'DYN_SS').gt.0)
 495       call reada(controlcard,"LIPSCALE",lipscale,1.0D0)
 496 allocate(ww(max_eneW))
 497       do i=1,n_eneW
 498         key = wname(i)(:ilen(wname(i)))
 499         call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
 500       enddo
 501
 502       write (iout,*) "iparm",iparm," myparm",myparm
 503 ! If reading not own parameters, skip assignment
 504
 505       if (iparm.eq.myparm .or. .not.separate_parset) then
 506
 507 !
 508 ! Setup weights for UNRES
 509 !
 510       wsc=ww(1)
 511       wscp=ww(2)
 512       welec=ww(3)
 513       wcorr=ww(4)
 514       wcorr5=ww(5)
 515       wcorr6=ww(6)
 516       wel_loc=ww(7)
 517       wturn3=ww(8)
 518       wturn4=ww(9)
 519       wturn6=ww(10)
 520       wang=ww(11)
 521       wscloc=ww(12)
 522       wtor=ww(13)
 523       wtor_d=ww(14)
 524       wstrain=ww(15)
 525       wvdwpp=ww(16)
 526       wbond=ww(18)
 527       wsccor=ww(19)
 528       wcatcat=ww(42)
 529       wcatprot=ww(41)
 530       wcorr3_nucl=ww(38)
 531       wcorr_nucl=ww(37)
 532       wtor_d_nucl=ww(36)
 533       wtor_nucl=ww(35)
 534       wsbloc=ww(34)
 535       wang_nucl=ww(33)
 536       wbond_nucl=ww(32)
 537       welsb=ww(31)
 538       wvdwsb=ww(30)
 539       welpsb=ww(29)
 540       wvdwpsb=ww(28)
 541       welpp=ww(27)
 542       wvdwpp_nucl=ww(26)
 543       wscbase=ww(46)
 544       wpepbase=ww(47)
 545       wscpho=ww(48)
 546       wpeppho=ww(49)
 547       wcatnucl=ww(50)
 548       wcat_tran=ww(56)
 549 !      print *,"KURWA",ww(48)
 550 !        "WSCBASE   ","WPEPBASE  ","WSCPHO    ","WPEPPHO   "
 551 !        "WVDWPP    ","WELPP     ","WVDWPSB   ","WELPSB    ","WVDWSB    ",&
 552 !        "WELSB     ","WBOND     ","WANG      ","WSBLOC    ","WTOR      ",&
 553 !        "WTORD     ","WCORR     ","WCORR3    ","WNULL     ","WNULL     ",&
 554 !        "WCATPROT  ","WCATCAT
 555 !       +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
 556 !       +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
 557 !       +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
 558 !       +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
 559
 560       endif
 561 !
 562 !el------
 563       allocate(weights(n_ene))
 564       weights(1)=wsc
 565       weights(2)=wscp
 566       weights(3)=welec
 567       weights(4)=wcorr
 568       weights(5)=wcorr5
 569       weights(6)=wcorr6
 570       weights(7)=wel_loc
 571       weights(8)=wturn3
 572       weights(9)=wturn4
 573       weights(10)=wturn6
 574       weights(11)=wang
 575       weights(12)=wscloc
 576       weights(13)=wtor
 577       weights(14)=wtor_d
 578       weights(15)=wstrain !0
 579       weights(16)=wvdwpp !
 580       weights(17)=wbond
 581       weights(18)=0 !scal14 !
 582       weights(21)=wsccor
 583       weights(42)=wcatprot
 584       weights(41)=wcatcat
 585       weights(26)=    wvdwpp_nucl
 586
 587       weights(27) =welpp
 588       weights(28) =wvdwpsb
 589       weights(29) =welpsb
 590       weights(30) =wvdwsb
 591       weights(31) =welsb
 592       weights(32) =wbond_nucl
 593       weights(33) =wang_nucl
 594       weights(34) =wsbloc
 595       weights(35) =wtor_nucl
 596       weights(36) =wtor_d_nucl
 597       weights(37) =wcorr_nucl
 598       weights(38) =wcorr3_nucl
 599       weights(41) =wcatcat
 600       weights(42) =wcatprot
 601       weights(46) =wscbase
 602       weights(47)= wpepbase
 603       weights(48) =wscpho
 604       weights(49) =wpeppho
 605       weights(50) =wcatnucl
 606       weights(56)=wcat_tran
 607
 608 ! el--------
 609       call card_concat(controlcard,.false.)
 610
 611 ! Return if not own parameters
 612
 613       if (iparm.ne.myparm .and. separate_parset) return
 614
 615       call reads(controlcard,"BONDPAR",bondname_t,bondname)
 616       open (ibond,file=bondname_t,status='old')
 617       rewind(ibond)
 618       call reads(controlcard,"THETPAR",thetname_t,thetname)
 619       open (ithep,file=thetname_t,status='old')
 620       rewind(ithep)
 621       call reads(controlcard,"ROTPAR",rotname_t,rotname)
 622       open (irotam,file=rotname_t,status='old')
 623       rewind(irotam)
 624       call reads(controlcard,"TORPAR",torname_t,torname)
 625       open (itorp,file=torname_t,status='old')
 626       rewind(itorp)
 627       call reads(controlcard,"TORDPAR",tordname_t,tordname)
 628       open (itordp,file=tordname_t,status='old')
 629       rewind(itordp)
 630       call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
 631       open (isccor,file=sccorname_t,status='old')
 632       rewind(isccor)
 633       call reads(controlcard,"FOURIER",fouriername_t,fouriername)
 634       open (ifourier,file=fouriername_t,status='old')
 635       rewind(ifourier)
 636       call reads(controlcard,"ELEPAR",elename_t,elename)
 637       open (ielep,file=elename_t,status='old')
 638       rewind(ielep)
 639       call reads(controlcard,"SIDEPAR",sidename_t,sidename)
 640       open (isidep,file=sidename_t,status='old')
 641       rewind(isidep)
 642       call reads(controlcard,"SCPPAR",scpname_t,scpname)
 643       open (iscpp,file=scpname_t,status='old')
 644       rewind(iscpp)
 645       call getenv_loc('IONPAR',ionname)
 646       open (iion,file=ionname,status='old')
 647       rewind(iion)
 648       write (iout,*) "Parameter set:",iparm
 649       write (iout,*) "Energy-term weights:"
 650       do i=1,n_eneW
 651         write (iout,'(i3,a16,f10.5)') i,wname(i),ww(i)
 652       enddo
 653       write (iout,*) "Sidechain potential file        : ",&
 654         sidename_t(:ilen(sidename_t))
 655 #ifndef OLDSCP
 656       write (iout,*) "SCp potential file              : ",&
 657         scpname_t(:ilen(scpname_t))
 658 #endif
 659       write (iout,*) "Electrostatic potential file    : ",&
 660         elename_t(:ilen(elename_t))
 661       write (iout,*) "Cumulant coefficient file       : ",&
 662         fouriername_t(:ilen(fouriername_t))
 663       write (iout,*) "Torsional parameter file        : ",&
 664         torname_t(:ilen(torname_t))
 665       write (iout,*) "Double torsional parameter file : ",&
 666         tordname_t(:ilen(tordname_t))
 667       write (iout,*) "Backbone-rotamer parameter file : ",&
 668         sccorname(:ilen(sccorname))
 669       write (iout,*) "Bond & inertia constant file    : ",&
 670         bondname_t(:ilen(bondname_t))
 671       write (iout,*) "Bending parameter file          : ",&
 672         thetname_t(:ilen(thetname_t))
 673       write (iout,*) "Rotamer parameter file          : ",&
 674         rotname_t(:ilen(rotname_t))
 675 #endif
 676 !
 677 ! Read the virtual-bond parameters, masses, and moments of inertia
 678 ! and Stokes' radii of the peptide group and side chains
 679 !
 680       allocate(dsc(ntyp1)) !(ntyp1)
 681       allocate(dsc_inv(ntyp1)) !(ntyp1)
 682       allocate(nbondterm(ntyp)) !(ntyp)
 683       allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 684       allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 685       allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
 686       allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
 687       allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
 688
 689 !el      allocate(msc(ntyp+1)) !(ntyp+1)
 690 !el      allocate(isc(ntyp+1)) !(ntyp+1)
 691 !el      allocate(restok(ntyp+1)) !(ntyp+1)
 692       allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 693
 694 #ifdef CRYST_BOND
 695       read (ibond,*) vbldp0,akp
 696       do i=1,ntyp
 697         nbondterm(i)=1
 698         read (ibond,*) vbldsc0(1,i),aksc(1,i)
 699         dsc(i) = vbldsc0(1,i)
 700         if (i.eq.10) then
 701           dsc_inv(i)=0.0D0
 702         else
 703           dsc_inv(i)=1.0D0/dsc(i)
 704         endif
 705       enddo
 706 #else
 707       read (ibond,*) ijunk,vbldp0,vbldpDUM,akp,rjunk
 708       do i=1,ntyp
 709         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
 710          j=1,nbondterm(i))
 711         dsc(i) = vbldsc0(1,i)
 712         if (i.eq.10) then
 713           dsc_inv(i)=0.0D0
 714         else
 715           dsc_inv(i)=1.0D0/dsc(i)
 716         endif
 717       enddo
 718 #endif
 719       if (lprint) then
 720         write(iout,'(/a/)')"Force constants virtual bonds:"
 721         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',&
 722          'inertia','Pstok'
 723         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
 724         do i=1,ntyp
 725           write (iout,'(a10,i3,6f10.5)') restyp(i,molnum(i)),nbondterm(i),&
 726             vbldsc0(1,i),aksc(1,i),abond0(1,i)
 727           do j=2,nbondterm(i)
 728             write (iout,'(13x,3f10.5)') &
 729               vbldsc0(j,i),aksc(j,i),abond0(j,i)
 730           enddo
 731         enddo
 732       endif
 733             if (.not. allocated(msc)) allocate(msc(-ntyp1:ntyp1,5))
 734             if (.not. allocated(restok)) allocate(restok(-ntyp1:ntyp1,5))
 735        if (oldion.eq.1) then
 736
 737             do i=1,ntyp_molec(5)
 738              read(iion,*) msc(i,5),restok(i,5)
 739              print *,msc(i,5),restok(i,5)
 740             enddo
 741             ip(5)=0.2
 742 !            isc(5)=0.2
 743             read (iion,*) ncatprotparm
 744             allocate(catprm(ncatprotparm,4))
 745             do k=1,4
 746             read (iion,*)  (catprm(i,k),i=1,ncatprotparm)
 747             enddo
 748             print *, catprm
 749       endif
 750       allocate(catnuclprm(14,ntyp_molec(2),ntyp_molec(5)))
 751       do i=1,ntyp_molec(5)
 752          do j=1,ntyp_molec(2)
 753          write(iout,*) i,j
 754             read(iionnucl,*) (catnuclprm(k,j,i),k=1,14)
 755          enddo
 756       enddo
 757       write(*,'(3(5x,a6)11(7x,a6))') "w1    ","w2    ","epslj ","pis1  ", &
 758       "sigma0","epsi0 ","chi1   ","chip1 ","sig   ","b1    ","b2    ", &
 759       "b3    ","b4    ","chis1  "
 760       do i=1,ntyp_molec(5)
 761          do j=1,ntyp_molec(2)
 762             write(*,'(3(f10.3,x),11(f12.6,x),a3,2a)') (catnuclprm(k,j,i),k=1,14), &
 763                                       restyp(i,5),"-",restyp(j,2)
 764          enddo
 765       enddo
 766
 767       read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
 768       do i=1,ntyp_molec(2)
 769         nbondterm_nucl(i)=1
 770         read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),rjunk,restok(i,2)
 771 !        dsc(i) = vbldsc0_nucl(1,i)
 772 !        if (i.eq.10) then
 773 !          dsc_inv(i)=0.0D0
 774 !        else
 775 !          dsc_inv(i)=1.0D0/dsc(i)
 776 !        endif
 777       enddo
 778
 779
 780 !----------------------------------------------------
 781       allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
 782       allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp))      !(-ntyp:ntyp)
 783       allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
 784       allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
 785       allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
 786       allocate(gthet(3,-ntyp:ntyp))     !(3,-ntyp:ntyp)
 787       do i=-ntyp,ntyp
 788         a0thet(i)=0.0D0
 789         do j=1,2
 790          do ichir1=-1,1
 791           do ichir2=-1,1
 792           athet(j,i,ichir1,ichir2)=0.0D0
 793           bthet(j,i,ichir1,ichir2)=0.0D0
 794           enddo
 795          enddo
 796         enddo
 797         do j=0,3
 798           polthet(j,i)=0.0D0
 799         enddo
 800         do j=1,3
 801           gthet(j,i)=0.0D0
 802         enddo
 803         theta0(i)=0.0D0
 804         sig0(i)=0.0D0
 805         sigc0(i)=0.0D0
 806       enddo
 807 !elwrite(iout,*) "parmread kontrol"
 808
 809 #ifdef CRYST_THETA
 810 !
 811 ! Read the parameters of the probability distribution/energy expression
 812 ! of the virtual-bond valence angles theta
 813 !
 814       do i=1,ntyp
 815         read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),&
 816           (bthet(j,i,1,1),j=1,2)
 817         read (ithep,*) (polthet(j,i),j=0,3)
 818 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 819         read (ithep,*) (gthet(j,i),j=1,3)
 820 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 821         read (ithep,*) theta0(i),sig0(i),sigc0(i)
 822         sigc0(i)=sigc0(i)**2
 823 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 824       enddo
 825 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 826       do i=1,ntyp
 827       athet(1,i,1,-1)=athet(1,i,1,1)
 828       athet(2,i,1,-1)=athet(2,i,1,1)
 829       bthet(1,i,1,-1)=-bthet(1,i,1,1)
 830       bthet(2,i,1,-1)=-bthet(2,i,1,1)
 831       athet(1,i,-1,1)=-athet(1,i,1,1)
 832       athet(2,i,-1,1)=-athet(2,i,1,1)
 833       bthet(1,i,-1,1)=bthet(1,i,1,1)
 834       bthet(2,i,-1,1)=bthet(2,i,1,1)
 835       enddo
 836 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 837       do i=-ntyp,-1
 838       a0thet(i)=a0thet(-i)
 839       athet(1,i,-1,-1)=athet(1,-i,1,1)
 840       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 841       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 842       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 843       athet(1,i,-1,1)=athet(1,-i,1,1)
 844       athet(2,i,-1,1)=-athet(2,-i,1,1)
 845       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 846       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 847       athet(1,i,1,-1)=-athet(1,-i,1,1)
 848       athet(2,i,1,-1)=athet(2,-i,1,1)
 849       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 850       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 851       theta0(i)=theta0(-i)
 852       sig0(i)=sig0(-i)
 853       sigc0(i)=sigc0(-i)
 854        do j=0,3
 855         polthet(j,i)=polthet(j,-i)
 856        enddo
 857        do j=1,3
 858          gthet(j,i)=gthet(j,-i)
 859        enddo
 860       enddo
 861 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 862       close (ithep)
 863 !elwrite(iout,*) "parmread kontrol in cryst_theta"
 864       if (lprint) then
 865 !       write (iout,'(a)')
 866 !    &    'Parameters of the virtual-bond valence angles:'
 867 !       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
 868 !    & '    ATHETA0   ','         A1   ','        A2    ',
 869 !    & '        B1    ','         B2   '
 870 !       do i=1,ntyp
 871 !         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 872 !    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
 873 !       enddo
 874 !       write (iout,'(/a/9x,5a/79(1h-))')
 875 !    & 'Parameters of the expression for sigma(theta_c):',
 876 !    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
 877 !    & '     ALPH3    ','    SIGMA0C   '
 878 !       do i=1,ntyp
 879 !         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
 880 !    &      (polthet(j,i),j=0,3),sigc0(i)
 881 !       enddo
 882 !       write (iout,'(/a/9x,5a/79(1h-))')
 883 !    & 'Parameters of the second gaussian:',
 884 !    & '    THETA0    ','     SIGMA0   ','        G1    ',
 885 !    & '        G2    ','         G3   '
 886 !       do i=1,ntyp
 887 !         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
 888 !    &       sig0(i),(gthet(j,i),j=1,3)
 889 !       enddo
 890         write (iout,'(a)') &
 891           'Parameters of the virtual-bond valence angles:'
 892         write (iout,'(/a/9x,5a/79(1h-))') &
 893        'Coefficients of expansion',&
 894        '     theta0   ','    a1*10^2   ','   a2*10^2    ',&
 895        '   b1*10^1    ','    b2*10^1   '
 896         do i=1,ntyp
 897           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
 898               a0thet(i),(100*athet(j,i,1,1),j=1,2),&
 899               (10*bthet(j,i,1,1),j=1,2)
 900         enddo
 901         write (iout,'(/a/9x,5a/79(1h-))') &
 902        'Parameters of the expression for sigma(theta_c):',&
 903        ' alpha0       ','  alph1       ',' alph2        ',&
 904        ' alhp3        ','   sigma0c    '
 905         do i=1,ntyp
 906           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
 907             (polthet(j,i),j=0,3),sigc0(i)
 908         enddo
 909         write (iout,'(/a/9x,5a/79(1h-))') &
 910        'Parameters of the second gaussian:',&
 911        '    theta0    ','  sigma0*10^2 ','      G1*10^-1',&
 912        '        G2    ','   G3*10^1    '
 913         do i=1,ntyp
 914           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
 915              100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 916         enddo
 917       endif
 918 #else
 919 !
 920 ! Read the parameters of Utheta determined from ab initio surfaces
 921 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 922 !
 923       allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
 924 !      write (iout,*) "tu dochodze"
 925       IF (tor_mode.eq.0) THEN
 926 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 927       read (ithep,*) nthetyp,ntheterm,ntheterm2,&
 928         ntheterm3,nsingle,ndouble
 929       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 930
 931 !----------------------------------------------------
 932 !      allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
 933       allocate(aa0thet(-nthetyp-1:nthetyp+1,&
 934         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 935 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 936       allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
 937         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 938 !(maxtheterm,-maxthetyp1:maxthetyp1,&
 939 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 940       allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 941         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 942       allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 943         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 944       allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 945         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 946       allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 947         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 948 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
 949 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 950       allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
 951         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 952       allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
 953         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 954 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
 955 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 956
 957
 958       read (ithep,*) (ithetyp(i),i=1,ntyp1)
 959       do i=-ntyp1,-1
 960         ithetyp(i)=-ithetyp(-i)
 961       enddo
 962 !      write (iout,*) "tu dochodze"
 963       aa0thet(:,:,:,:)=0.0d0
 964       aathet(:,:,:,:,:)=0.0d0
 965       bbthet(:,:,:,:,:,:)=0.0d0
 966       ccthet(:,:,:,:,:,:)=0.0d0
 967       ddthet(:,:,:,:,:,:)=0.0d0
 968       eethet(:,:,:,:,:,:)=0.0d0
 969       ffthet(:,:,:,:,:,:,:)=0.0d0
 970       ggthet(:,:,:,:,:,:,:)=0.0d0
 971
 972       do iblock=1,2
 973       do i=0,nthetyp
 974         do j=-nthetyp,nthetyp
 975           do k=-nthetyp,nthetyp
 976             read (ithep,'(6a)') res1
 977             read (ithep,*) aa0thet(i,j,k,iblock)
 978             read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
 979             read (ithep,*) &
 980              ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
 981               (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
 982               (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
 983               (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
 984               ll=1,ntheterm2)
 985             read (ithep,*) &
 986             (((ffthet(llll,lll,ll,i,j,k,iblock),&
 987                ffthet(lll,llll,ll,i,j,k,iblock),&
 988                ggthet(llll,lll,ll,i,j,k,iblock),&
 989                ggthet(lll,llll,ll,i,j,k,iblock),&
 990                llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
 991           enddo
 992         enddo
 993       enddo
 994 !
 995 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
 996 ! coefficients of theta-and-gamma-dependent terms are zero.
 997 !
 998       do i=1,nthetyp
 999         do j=1,nthetyp
1000           do l=1,ntheterm
1001             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1002             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1003           enddo
1004           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1005           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1006         enddo
1007         do l=1,ntheterm
1008           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1009         enddo
1010         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1011       enddo
1012       enddo
1013 ! Substitution for D aminoacids from symmetry.
1014       do iblock=1,2
1015       do i=-nthetyp,0
1016         do j=-nthetyp,nthetyp
1017           do k=-nthetyp,nthetyp
1018            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1019            do l=1,ntheterm
1020            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1021            enddo
1022            do ll=1,ntheterm2
1023             do lll=1,nsingle
1024             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1025             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1026             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1027             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1028             enddo
1029           enddo
1030           do ll=1,ntheterm3
1031            do lll=2,ndouble
1032             do llll=1,lll-1
1033             ffthet(llll,lll,ll,i,j,k,iblock)= &
1034             ffthet(llll,lll,ll,-i,-j,-k,iblock)
1035             ffthet(lll,llll,ll,i,j,k,iblock)= &
1036             ffthet(lll,llll,ll,-i,-j,-k,iblock)
1037             ggthet(llll,lll,ll,i,j,k,iblock)= &
1038             -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1039             ggthet(lll,llll,ll,i,j,k,iblock)= &
1040             -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1041             enddo !ll
1042            enddo  !lll
1043           enddo   !llll
1044          enddo    !k
1045         enddo     !j
1046        enddo      !i
1047       enddo       !iblock
1048
1049 !
1050 ! Control printout of the coefficients of virtual-bond-angle potentials
1051 !
1052       do iblock=1,2
1053       if (lprint) then
1054         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1055         do i=1,nthetyp+1
1056           do j=1,nthetyp+1
1057             do k=1,nthetyp+1
1058               write (iout,'(//4a)') &
1059                'Type ',onelett(i),onelett(j),onelett(k)
1060               write (iout,'(//a,10x,a)') " l","a[l]"
1061               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1062               write (iout,'(i2,1pe15.5)') &
1063                  (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1064             do l=1,ntheterm2
1065               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))') &
1066                 "b",l,"c",l,"d",l,"e",l
1067               do m=1,nsingle
1068                 write (iout,'(i2,4(1pe15.5))') m,&
1069                 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1070                 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1071               enddo
1072             enddo
1073             do l=1,ntheterm3
1074               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))') &
1075                 "f+",l,"f-",l,"g+",l,"g-",l
1076               do m=2,ndouble
1077                 do n=1,m-1
1078                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1079                     ffthet(n,m,l,i,j,k,iblock),&
1080                     ffthet(m,n,l,i,j,k,iblock),&
1081                     ggthet(n,m,l,i,j,k,iblock),&
1082                     ggthet(m,n,l,i,j,k,iblock)
1083                 enddo
1084               enddo
1085             enddo
1086           enddo
1087         enddo
1088       enddo
1089       call flush(iout)
1090       endif
1091       enddo
1092       ELSE
1093 !C here will be the apropriate recalibrating for D-aminoacid
1094       read (ithep,*) nthetyp
1095       allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1096       allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1097       do i=-nthetyp+1,nthetyp-1
1098         read (ithep,*) nbend_kcc_Tb(i)
1099         do j=0,nbend_kcc_Tb(i)
1100           read (ithep,*) ijunk,v1bend_chyb(j,i)
1101         enddo
1102       enddo
1103       if (lprint) then
1104         write (iout,'(a)') &
1105          "Parameters of the valence-only potentials"
1106         do i=-nthetyp+1,nthetyp-1
1107         write (iout,'(2a)') "Type ",toronelet(i)
1108         do k=0,nbend_kcc_Tb(i)
1109           write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1110         enddo
1111         enddo
1112       endif
1113       ENDIF ! TOR_MODE
1114
1115 #endif
1116 !--------------- Reading theta parameters for nucleic acid-------
1117       read (ithep_nucl,*) nthetyp_nucl,ntheterm_nucl,&
1118       ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1119       nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1120       allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1121       allocate(aa0thet_nucl(nthetyp_nucl+1,&
1122         nthetyp_nucl+1,nthetyp_nucl+1))
1123 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1124       allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1125         nthetyp_nucl+1,nthetyp_nucl+1))
1126 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1127 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1128       allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1129         nthetyp_nucl+1,nthetyp_nucl+1))
1130       allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1131         nthetyp_nucl+1,nthetyp_nucl+1))
1132       allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1133         nthetyp_nucl+1,nthetyp_nucl+1))
1134       allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1135         nthetyp_nucl+1,nthetyp_nucl+1))
1136 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1137 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1138       allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1139          nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1140       allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1141          nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1142
1143 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1144 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1145
1146       read (ithep_nucl,*) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1147
1148       aa0thet_nucl(:,:,:)=0.0d0
1149       aathet_nucl(:,:,:,:)=0.0d0
1150       bbthet_nucl(:,:,:,:,:)=0.0d0
1151       ccthet_nucl(:,:,:,:,:)=0.0d0
1152       ddthet_nucl(:,:,:,:,:)=0.0d0
1153       eethet_nucl(:,:,:,:,:)=0.0d0
1154       ffthet_nucl(:,:,:,:,:,:)=0.0d0
1155       ggthet_nucl(:,:,:,:,:,:)=0.0d0
1156
1157       do i=1,nthetyp_nucl
1158         do j=1,nthetyp_nucl
1159           do k=1,nthetyp_nucl
1160             read (ithep_nucl,'(3a)') t1,t2,t3
1161             read (ithep_nucl,*) aa0thet_nucl(i,j,k)
1162             read (ithep_nucl,*)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1163             read (ithep_nucl,*) &
1164             (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1165             (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1166             (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1167             (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1168             read (ithep_nucl,*) &
1169            (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1170               ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1171               llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1172           enddo
1173         enddo
1174       enddo
1175
1176
1177 !-------------------------------------------
1178       allocate(nlob(ntyp1)) !(ntyp1)
1179       allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1180       allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1181       allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1182
1183       do i=1,ntyp
1184         do j=1,maxlob
1185           bsc(j,i)=0.0D0
1186           nlob(i)=0
1187         enddo
1188       enddo
1189       nlob(ntyp1)=0
1190       dsc(ntyp1)=0.0D0
1191
1192       do i=-ntyp,ntyp
1193         do j=1,maxlob
1194           do k=1,3
1195             censc(k,j,i)=0.0D0
1196           enddo
1197           do k=1,3
1198             do l=1,3
1199               gaussc(l,k,j,i)=0.0D0
1200             enddo
1201           enddo
1202         enddo
1203       enddo
1204
1205 #ifdef CRYST_SC
1206 !
1207 ! Read the parameters of the probability distribution/energy expression
1208 ! of the side chains.
1209 !
1210       do i=1,ntyp
1211 !c      write (iout,*) "tu dochodze",i
1212         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
1213         if (i.eq.10) then
1214           dsc_inv(i)=0.0D0
1215         else
1216           dsc_inv(i)=1.0D0/dsc(i)
1217         endif
1218         if (i.ne.10) then
1219         do j=1,nlob(i)
1220           do k=1,3
1221             do l=1,3
1222               blower(l,k,j)=0.0D0
1223             enddo
1224           enddo
1225         enddo
1226         bsc(1,i)=0.0D0
1227         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
1228         censc(1,1,-i)=censc(1,1,i)
1229         censc(2,1,-i)=censc(2,1,i)
1230         censc(3,1,-i)=-censc(3,1,i)
1231         do j=2,nlob(i)
1232           read (irotam,*) bsc(j,i)
1233           read (irotam,*) (censc(k,j,i),k=1,3),&
1234                                        ((blower(k,l,j),l=1,k),k=1,3)
1235         censc(1,j,-i)=censc(1,j,i)
1236         censc(2,j,-i)=censc(2,j,i)
1237         censc(3,j,-i)=-censc(3,j,i)
1238 ! BSC is amplitude of Gaussian
1239         enddo
1240         do j=1,nlob(i)
1241           do k=1,3
1242             do l=1,k
1243               akl=0.0D0
1244               do m=1,3
1245                 akl=akl+blower(k,m,j)*blower(l,m,j)
1246               enddo
1247               gaussc(k,l,j,i)=akl
1248               gaussc(l,k,j,i)=akl
1249              if (((k.eq.3).and.(l.ne.3)) &
1250               .or.((l.eq.3).and.(k.ne.3))) then
1251                 gaussc(k,l,j,-i)=-akl
1252                 gaussc(l,k,j,-i)=-akl
1253               else
1254                 gaussc(k,l,j,-i)=akl
1255                 gaussc(l,k,j,-i)=akl
1256               endif
1257             enddo
1258           enddo
1259         enddo
1260         endif
1261       enddo
1262       close (irotam)
1263       if (lprint) then
1264         write (iout,'(/a)') 'Parameters of side-chain local geometry'
1265         do i=1,ntyp
1266           nlobi=nlob(i)
1267           if (nlobi.gt.0) then
1268           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
1269            ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1270 !          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1271 !          write (iout,'(a,f10.4,4(16x,f10.4))')
1272 !     &                             'Center  ',(bsc(j,i),j=1,nlobi)
1273 !          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
1274            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1275                                    'log h',(bsc(j,i),j=1,nlobi)
1276            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1277           'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1278 !          write (iout,'(a)')
1279 !         do j=1,nlobi
1280 !           ind=0
1281 !           do k=1,3
1282 !             do l=1,k
1283 !              ind=ind+1
1284 !              blower(k,l,j)=gaussc(ind,j,i)
1285 !             enddo
1286 !           enddo
1287 !         enddo
1288           do k=1,3
1289             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1290                        ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1291           enddo
1292           endif
1293         enddo
1294       endif
1295 #else
1296 !
1297 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1298 ! added by Urszula Kozlowska 07/11/2007
1299 !
1300       allocate(sc_parmin(65,ntyp))      !(maxsccoef,ntyp)
1301
1302       do i=1,ntyp
1303         read (irotam,*)
1304        if (i.eq.10) then
1305          read (irotam,*)
1306        else
1307          do j=1,65
1308            read(irotam,*) sc_parmin(j,i)
1309          enddo
1310        endif
1311       enddo
1312 #endif
1313       close(irotam)
1314 !---------reading nucleic acid parameters for rotamers-------------------
1315       allocate(sc_parmin_nucl(9,ntyp_molec(2)))      !(maxsccoef,ntyp)
1316       do i=1,ntyp_molec(2)
1317         read (irotam_nucl,*)
1318         do j=1,9
1319           read(irotam_nucl,*) sc_parmin_nucl(j,i)
1320         enddo
1321       enddo
1322       close(irotam_nucl)
1323       if (lprint) then
1324         write (iout,*)
1325         write (iout,*) "Base rotamer parameters"
1326         do i=1,ntyp_molec(2)
1327           write (iout,'(a)') restyp(i,2)
1328           write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1329         enddo
1330       endif
1331
1332
1333       read (ifourier,*) nloctyp
1334 !el write(iout,*)"nloctyp",nloctyp
1335       SPLIT_FOURIERTOR = nloctyp.lt.0
1336       nloctyp = iabs(nloctyp)
1337 #ifdef NEWCORR
1338       if (.not.allocated(itype2loc)) allocate(itype2loc(-ntyp1:ntyp1))
1339        print *,"shape",shape(itype2loc)
1340       allocate(iloctyp(-nloctyp:nloctyp))
1341       allocate(bnew1(3,2,-nloctyp:nloctyp))
1342       allocate(bnew2(3,2,-nloctyp:nloctyp))
1343       allocate(ccnew(3,2,-nloctyp:nloctyp))
1344       allocate(ddnew(3,2,-nloctyp:nloctyp))
1345       allocate(e0new(3,-nloctyp:nloctyp))
1346       allocate(eenew(2,2,2,-nloctyp:nloctyp))
1347       allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1348       allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1349       allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1350       allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1351       allocate(e0newtor(3,-nloctyp:nloctyp))
1352       allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1353
1354       read (ifourier,*) (itype2loc(i),i=1,ntyp)
1355       read (ifourier,*) (iloctyp(i),i=0,nloctyp-1)
1356       itype2loc(ntyp1)=nloctyp
1357       iloctyp(nloctyp)=ntyp1
1358       do i=1,ntyp1
1359         itype2loc(-i)=-itype2loc(i)
1360       enddo
1361 #else
1362       allocate(iloctyp(-nloctyp:nloctyp))
1363       allocate(itype2loc(-ntyp1:ntyp1))
1364       iloctyp(0)=10
1365       iloctyp(1)=9
1366       iloctyp(2)=20
1367       iloctyp(3)=ntyp1
1368 #endif
1369       do i=1,nloctyp
1370         iloctyp(-i)=-iloctyp(i)
1371       enddo
1372 !c      write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1373 !c      write (iout,*) "nloctyp",nloctyp,
1374 !c     &  " iloctyp",(iloctyp(i),i=0,nloctyp)
1375 !c      write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1376 !c      write (iout,*) "nloctyp",nloctyp,
1377 !c     &  " iloctyp",(iloctyp(i),i=0,nloctyp)
1378 #ifdef NEWCORR
1379       do i=0,nloctyp-1
1380 !c             write (iout,*) "NEWCORR",i
1381         read (ifourier,*)
1382         do ii=1,3
1383           do j=1,2
1384             read (ifourier,*) bnew1(ii,j,i)
1385           enddo
1386         enddo
1387 !c             write (iout,*) "NEWCORR BNEW1"
1388 !c             write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1389         do ii=1,3
1390           do j=1,2
1391             read (ifourier,*) bnew2(ii,j,i)
1392           enddo
1393         enddo
1394 !c             write (iout,*) "NEWCORR BNEW2"
1395 !c             write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1396         do kk=1,3
1397           read (ifourier,*) ccnew(kk,1,i)
1398           read (ifourier,*) ccnew(kk,2,i)
1399         enddo
1400 !c             write (iout,*) "NEWCORR CCNEW"
1401 !c             write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1402         do kk=1,3
1403           read (ifourier,*) ddnew(kk,1,i)
1404           read (ifourier,*) ddnew(kk,2,i)
1405         enddo
1406 !c             write (iout,*) "NEWCORR DDNEW"
1407 !c             write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1408         do ii=1,2
1409           do jj=1,2
1410             do kk=1,2
1411               read (ifourier,*) eenew(ii,jj,kk,i)
1412             enddo
1413           enddo
1414         enddo
1415 !c             write (iout,*) "NEWCORR EENEW1"
1416 !c             write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1417         do ii=1,3
1418           read (ifourier,*) e0new(ii,i)
1419         enddo
1420 !c             write (iout,*) (e0new(ii,i),ii=1,3)
1421       enddo
1422 !c             write (iout,*) "NEWCORR EENEW"
1423       do i=0,nloctyp-1
1424         do ii=1,3
1425           ccnew(ii,1,i)=ccnew(ii,1,i)/2
1426           ccnew(ii,2,i)=ccnew(ii,2,i)/2
1427           ddnew(ii,1,i)=ddnew(ii,1,i)/2
1428           ddnew(ii,2,i)=ddnew(ii,2,i)/2
1429         enddo
1430       enddo
1431       do i=1,nloctyp-1
1432         do ii=1,3
1433           bnew1(ii,1,-i)= bnew1(ii,1,i)
1434           bnew1(ii,2,-i)=-bnew1(ii,2,i)
1435           bnew2(ii,1,-i)= bnew2(ii,1,i)
1436           bnew2(ii,2,-i)=-bnew2(ii,2,i)
1437         enddo
1438         do ii=1,3
1439 !c          ccnew(ii,1,i)=ccnew(ii,1,i)/2
1440 !c          ccnew(ii,2,i)=ccnew(ii,2,i)/2
1441 !c          ddnew(ii,1,i)=ddnew(ii,1,i)/2
1442 !c          ddnew(ii,2,i)=ddnew(ii,2,i)/2
1443           ccnew(ii,1,-i)=ccnew(ii,1,i)
1444           ccnew(ii,2,-i)=-ccnew(ii,2,i)
1445           ddnew(ii,1,-i)=ddnew(ii,1,i)
1446           ddnew(ii,2,-i)=-ddnew(ii,2,i)
1447         enddo
1448         e0new(1,-i)= e0new(1,i)
1449         e0new(2,-i)=-e0new(2,i)
1450         e0new(3,-i)=-e0new(3,i)
1451         do kk=1,2
1452           eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1453           eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1454           eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1455           eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1456         enddo
1457       enddo
1458       if (lprint) then
1459         write (iout,'(a)') "Coefficients of the multibody terms"
1460         do i=-nloctyp+1,nloctyp-1
1461           write (iout,*) "Type: ",onelet(iloctyp(i))
1462           write (iout,*) "Coefficients of the expansion of B1"
1463           do j=1,2
1464             write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1465           enddo
1466           write (iout,*) "Coefficients of the expansion of B2"
1467           do j=1,2
1468             write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1469           enddo
1470           write (iout,*) "Coefficients of the expansion of C"
1471           write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1472           write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1473           write (iout,*) "Coefficients of the expansion of D"
1474           write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1475           write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1476           write (iout,*) "Coefficients of the expansion of E"
1477           write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1478           do j=1,2
1479             do k=1,2
1480               write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1481             enddo
1482           enddo
1483         enddo
1484       endif
1485       IF (SPLIT_FOURIERTOR) THEN
1486       do i=0,nloctyp-1
1487 !c             write (iout,*) "NEWCORR TOR",i
1488         read (ifourier,*)
1489         do ii=1,3
1490           do j=1,2
1491             read (ifourier,*) bnew1tor(ii,j,i)
1492           enddo
1493         enddo
1494 !c             write (iout,*) "NEWCORR BNEW1 TOR"
1495 !c             write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1496         do ii=1,3
1497           do j=1,2
1498             read (ifourier,*) bnew2tor(ii,j,i)
1499           enddo
1500         enddo
1501 !c             write (iout,*) "NEWCORR BNEW2 TOR"
1502 !c             write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1503         do kk=1,3
1504           read (ifourier,*) ccnewtor(kk,1,i)
1505           read (ifourier,*) ccnewtor(kk,2,i)
1506         enddo
1507 !c             write (iout,*) "NEWCORR CCNEW TOR"
1508 !c             write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1509         do kk=1,3
1510           read (ifourier,*) ddnewtor(kk,1,i)
1511           read (ifourier,*) ddnewtor(kk,2,i)
1512         enddo
1513 !c             write (iout,*) "NEWCORR DDNEW TOR"
1514 !c             write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1515         do ii=1,2
1516           do jj=1,2
1517             do kk=1,2
1518               read (ifourier,*) eenewtor(ii,jj,kk,i)
1519             enddo
1520           enddo
1521         enddo
1522 !c         write (iout,*) "NEWCORR EENEW1 TOR"
1523 !c         write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1524         do ii=1,3
1525           read (ifourier,*) e0newtor(ii,i)
1526         enddo
1527 !c             write (iout,*) (e0newtor(ii,i),ii=1,3)
1528       enddo
1529 !c             write (iout,*) "NEWCORR EENEW TOR"
1530       do i=0,nloctyp-1
1531         do ii=1,3
1532           ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1533           ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1534           ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1535           ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1536         enddo
1537       enddo
1538       do i=1,nloctyp-1
1539         do ii=1,3
1540           bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1541           bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1542           bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1543           bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1544         enddo
1545         do ii=1,3
1546           ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1547           ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1548           ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1549           ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1550         enddo
1551         e0newtor(1,-i)= e0newtor(1,i)
1552         e0newtor(2,-i)=-e0newtor(2,i)
1553         e0newtor(3,-i)=-e0newtor(3,i)
1554         do kk=1,2
1555           eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1556           eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1557           eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1558           eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1559         enddo
1560       enddo
1561       if (lprint) then
1562         write (iout,'(a)') &
1563          "Single-body coefficients of the torsional potentials"
1564         do i=-nloctyp+1,nloctyp-1
1565           write (iout,*) "Type: ",onelet(iloctyp(i))
1566           write (iout,*) "Coefficients of the expansion of B1tor"
1567           do j=1,2
1568             write (iout,'(3hB1(,i1,1h),3f10.5)') &
1569              j,(bnew1tor(k,j,i),k=1,3)
1570           enddo
1571           write (iout,*) "Coefficients of the expansion of B2tor"
1572           do j=1,2
1573             write (iout,'(3hB2(,i1,1h),3f10.5)') &
1574              j,(bnew2tor(k,j,i),k=1,3)
1575           enddo
1576           write (iout,*) "Coefficients of the expansion of Ctor"
1577           write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1578           write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1579           write (iout,*) "Coefficients of the expansion of Dtor"
1580           write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1581           write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1582           write (iout,*) "Coefficients of the expansion of Etor"
1583           write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1584           do j=1,2
1585             do k=1,2
1586               write (iout,'(1hE,2i1,2f10.5)') &
1587                j,k,(eenewtor(l,j,k,i),l=1,2)
1588             enddo
1589           enddo
1590         enddo
1591       endif
1592       ELSE
1593       do i=-nloctyp+1,nloctyp-1
1594         do ii=1,3
1595           do j=1,2
1596             bnew1tor(ii,j,i)=bnew1(ii,j,i)
1597           enddo
1598         enddo
1599         do ii=1,3
1600           do j=1,2
1601             bnew2tor(ii,j,i)=bnew2(ii,j,i)
1602           enddo
1603         enddo
1604         do ii=1,3
1605           ccnewtor(ii,1,i)=ccnew(ii,1,i)
1606           ccnewtor(ii,2,i)=ccnew(ii,2,i)
1607           ddnewtor(ii,1,i)=ddnew(ii,1,i)
1608           ddnewtor(ii,2,i)=ddnew(ii,2,i)
1609         enddo
1610       enddo
1611       ENDIF !SPLIT_FOURIER_TOR
1612 #else
1613       allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1614       allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1615       allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1616       allocate(b(13,-nloctyp-1:nloctyp+1))
1617       if (lprint) &
1618        write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1619       do i=0,nloctyp-1
1620         read (ifourier,*)
1621         read (ifourier,*) (b(ii,i),ii=1,13)
1622         if (lprint) then
1623         write (iout,*) 'Type ',onelet(iloctyp(i))
1624         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1625         endif
1626         if (i.gt.0) then
1627         b(2,-i)= b(2,i)
1628         b(3,-i)= b(3,i)
1629         b(4,-i)=-b(4,i)
1630         b(5,-i)=-b(5,i)
1631         endif
1632         CCold(1,1,i)= b(7,i)
1633         CCold(2,2,i)=-b(7,i)
1634         CCold(2,1,i)= b(9,i)
1635         CCold(1,2,i)= b(9,i)
1636         CCold(1,1,-i)= b(7,i)
1637         CCold(2,2,-i)=-b(7,i)
1638         CCold(2,1,-i)=-b(9,i)
1639         CCold(1,2,-i)=-b(9,i)
1640         DDold(1,1,i)= b(6,i)
1641         DDold(2,2,i)=-b(6,i)
1642         DDold(2,1,i)= b(8,i)
1643         DDold(1,2,i)= b(8,i)
1644         DDold(1,1,-i)= b(6,i)
1645         DDold(2,2,-i)=-b(6,i)
1646         DDold(2,1,-i)=-b(8,i)
1647         DDold(1,2,-i)=-b(8,i)
1648         EEold(1,1,i)= b(10,i)+b(11,i)
1649         EEold(2,2,i)=-b(10,i)+b(11,i)
1650         EEold(2,1,i)= b(12,i)-b(13,i)
1651         EEold(1,2,i)= b(12,i)+b(13,i)
1652         EEold(1,1,-i)= b(10,i)+b(11,i)
1653         EEold(2,2,-i)=-b(10,i)+b(11,i)
1654         EEold(2,1,-i)=-b(12,i)+b(13,i)
1655         EEold(1,2,-i)=-b(12,i)-b(13,i)
1656         write(iout,*) "TU DOCHODZE"
1657         print *,"JESTEM"
1658       enddo
1659       if (lprint) then
1660       write (iout,*)
1661       write (iout,*) &
1662       "Coefficients of the cumulants (independent of valence angles)"
1663       do i=-nloctyp+1,nloctyp-1
1664         write (iout,*) 'Type ',onelet(iloctyp(i))
1665         write (iout,*) 'B1'
1666         write(iout,'(2f10.5)') B(3,i),B(5,i)
1667         write (iout,*) 'B2'
1668         write(iout,'(2f10.5)') B(2,i),B(4,i)
1669         write (iout,*) 'CC'
1670         do j=1,2
1671           write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1672         enddo
1673         write(iout,*) 'DD'
1674         do j=1,2
1675           write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1676         enddo
1677         write(iout,*) 'EE'
1678         do j=1,2
1679           write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1680         enddo
1681       enddo
1682       endif
1683 #endif
1684 #ifdef CRYST_TOR
1685 !
1686 ! Read torsional parameters in old format
1687 !
1688       allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1689
1690       read (itorp,*) ntortyp,nterm_old
1691       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1692       read (itorp,*) (itortyp(i),i=1,ntyp)
1693
1694 !el from energy module--------
1695       allocate(v1(nterm_old,ntortyp,ntortyp))
1696       allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1697 !el---------------------------
1698
1699       do i=1,ntortyp
1700         do j=1,ntortyp
1701           read (itorp,'(a)')
1702           do k=1,nterm_old
1703             read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
1704           enddo
1705         enddo
1706       enddo
1707       close (itorp)
1708       if (lprint) then
1709         write (iout,'(/a/)') 'Torsional constants:'
1710         do i=1,ntortyp
1711           do j=1,ntortyp
1712             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1713             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1714           enddo
1715         enddo
1716       endif
1717
1718
1719 #else
1720 !
1721 ! Read torsional parameters
1722 !
1723       IF (TOR_MODE.eq.0) THEN
1724
1725       allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1726
1727       read (itorp,*) ntortyp
1728       read (itorp,*) (itortyp(i),i=1,ntyp)
1729       write (iout,*) 'ntortyp',ntortyp
1730
1731 !el from energy module---------
1732       allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1733       allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1734
1735       allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1736       allocate(vlor2(maxlor,ntortyp,ntortyp))
1737       allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1738       allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1739
1740       allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1741       allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1742 !el---------------------------
1743       do iblock=1,2
1744         do i=-ntortyp,ntortyp
1745           do j=-ntortyp,ntortyp
1746             nterm(i,j,iblock)=0
1747             nlor(i,j,iblock)=0
1748           enddo
1749         enddo
1750       enddo
1751 !el---------------------------
1752
1753       do iblock=1,2
1754       do i=-ntyp,-1
1755        itortyp(i)=-itortyp(-i)
1756       enddo
1757 !      write (iout,*) 'ntortyp',ntortyp
1758       do i=0,ntortyp-1
1759         do j=-ntortyp+1,ntortyp-1
1760           read (itorp,*) nterm(i,j,iblock),&
1761                 nlor(i,j,iblock)
1762           nterm(-i,-j,iblock)=nterm(i,j,iblock)
1763           nlor(-i,-j,iblock)=nlor(i,j,iblock)
1764           v0ij=0.0d0
1765           si=-1.0d0
1766           do k=1,nterm(i,j,iblock)
1767             read (itorp,*) kk,v1(k,i,j,iblock),&
1768             v2(k,i,j,iblock)
1769             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1770             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1771             v0ij=v0ij+si*v1(k,i,j,iblock)
1772             si=-si
1773          enddo
1774           do k=1,nlor(i,j,iblock)
1775             read (itorp,*) kk,vlor1(k,i,j),&
1776               vlor2(k,i,j),vlor3(k,i,j)
1777             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1778           enddo
1779           v0(i,j,iblock)=v0ij
1780           v0(-i,-j,iblock)=v0ij
1781         enddo
1782       enddo
1783       enddo
1784       close (itorp)
1785       if (lprint) then
1786         do iblock=1,2 !el
1787         write (iout,'(/a/)') 'Torsional constants:'
1788         do i=1,ntortyp
1789           do j=1,ntortyp
1790             write (iout,*) 'ityp',i,' jtyp',j
1791             write (iout,*) 'Fourier constants'
1792             do k=1,nterm(i,j,iblock)
1793               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1794               v2(k,i,j,iblock)
1795             enddo
1796             write (iout,*) 'Lorenz constants'
1797             do k=1,nlor(i,j,iblock)
1798               write (iout,'(3(1pe15.5))') &
1799                vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1800             enddo
1801           enddo
1802         enddo
1803         enddo
1804       endif
1805 !
1806 ! 6/23/01 Read parameters for double torsionals
1807 !
1808 !el from energy module------------
1809       allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1810       allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1811 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1812       allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1813       allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1814         !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1815       allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1816       allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1817         !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1818 !---------------------------------
1819
1820       do iblock=1,2
1821       do i=0,ntortyp-1
1822         do j=-ntortyp+1,ntortyp-1
1823           do k=-ntortyp+1,ntortyp-1
1824             read (itordp,'(3a1)') t1,t2,t3
1825 !              write (iout,*) "OK onelett",
1826 !     &         i,j,k,t1,t2,t3
1827
1828             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1829               .or. t3.ne.toronelet(k)) then
1830               write (iout,*) "Error in double torsional parameter file",&
1831                i,j,k,t1,t2,t3
1832 #ifdef MPI
1833               call MPI_Finalize(Ierror)
1834 #endif
1835                stop "Error in double torsional parameter file"
1836             endif
1837           read (itordp,*) ntermd_1(i,j,k,iblock),&
1838                ntermd_2(i,j,k,iblock)
1839             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1840             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1841             read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,&
1842                ntermd_1(i,j,k,iblock))
1843             read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,&
1844                ntermd_1(i,j,k,iblock))
1845             read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,&
1846                ntermd_1(i,j,k,iblock))
1847             read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,&
1848                ntermd_1(i,j,k,iblock))
1849 ! Martix of D parameters for one dimesional foureir series
1850             do l=1,ntermd_1(i,j,k,iblock)
1851              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1852              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1853              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1854              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1855 !            write(iout,*) "whcodze" ,
1856 !     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1857             enddo
1858             read (itordp,*) ((v2c(l,m,i,j,k,iblock),&
1859                v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1860                v2s(m,l,i,j,k,iblock),&
1861                m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1862 ! Martix of D parameters for two dimesional fourier series
1863             do l=1,ntermd_2(i,j,k,iblock)
1864              do m=1,l-1
1865              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1866              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1867              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1868              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1869              enddo!m
1870             enddo!l
1871           enddo!k
1872         enddo!j
1873       enddo!i
1874       enddo!iblock
1875       if (lprint) then
1876       write (iout,*)
1877       write (iout,*) 'Constants for double torsionals'
1878       do iblock=1,2
1879       do i=0,ntortyp-1
1880         do j=-ntortyp+1,ntortyp-1
1881           do k=-ntortyp+1,ntortyp-1
1882             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1883               ' nsingle',ntermd_1(i,j,k,iblock),&
1884               ' ndouble',ntermd_2(i,j,k,iblock)
1885             write (iout,*)
1886             write (iout,*) 'Single angles:'
1887             do l=1,ntermd_1(i,j,k,iblock)
1888               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1889                  v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1890                  v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1891                  v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1892             enddo
1893             write (iout,*)
1894             write (iout,*) 'Pairs of angles:'
1895             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1896             do l=1,ntermd_2(i,j,k,iblock)
1897               write (iout,'(i5,20f10.5)') &
1898                l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1899             enddo
1900             write (iout,*)
1901            write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1902             do l=1,ntermd_2(i,j,k,iblock)
1903               write (iout,'(i5,20f10.5)') &
1904                l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1905                (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1906             enddo
1907             write (iout,*)
1908           enddo
1909         enddo
1910       enddo
1911       enddo
1912       endif
1913 #ifndef NEWCORR
1914       do i=1,ntyp1
1915         itype2loc(i)=itortyp(i)
1916       enddo
1917 #endif
1918       ELSE IF (TOR_MODE.eq.1) THEN
1919
1920 !C read valence-torsional parameters
1921       read (itorp,*) ntortyp
1922       nkcctyp=ntortyp
1923       write (iout,*) "Valence-torsional parameters read in ntortyp",&
1924         ntortyp
1925       read (itorp,*) (itortyp(i),i=1,ntyp)
1926       write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
1927 #ifndef NEWCORR
1928       do i=1,ntyp1
1929         itype2loc(i)=itortyp(i)
1930       enddo
1931 #endif
1932       do i=-ntyp,-1
1933         itortyp(i)=-itortyp(-i)
1934       enddo
1935       do i=-ntortyp+1,ntortyp-1
1936         do j=-ntortyp+1,ntortyp-1
1937 !C first we read the cos and sin gamma parameters
1938           read (itorp,'(13x,a)') string
1939           write (iout,*) i,j,string
1940           read (itorp,*) &
1941          nterm_kcc(j,i),nterm_kcc_Tb(j,i)
1942 !C           read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
1943           do k=1,nterm_kcc(j,i)
1944             do l=1,nterm_kcc_Tb(j,i)
1945               do ll=1,nterm_kcc_Tb(j,i)
1946               read (itorp,*) ii,jj,kk, &
1947                v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
1948               enddo
1949             enddo
1950           enddo
1951         enddo
1952       enddo
1953       ELSE
1954 #ifdef NEWCORR
1955 !c AL 4/8/16: Calculate coefficients from one-body parameters
1956       ntortyp=nloctyp
1957       allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1958       allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
1959       allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
1960       allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
1961       allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
1962
1963       do i=-ntyp1,ntyp1
1964        print *,i,itortyp(i)
1965        itortyp(i)=itype2loc(i)
1966       enddo
1967       write (iout,*) &
1968       "Val-tor parameters calculated from cumulant coefficients ntortyp"&
1969       ,ntortyp
1970       do i=-ntortyp+1,ntortyp-1
1971         do j=-ntortyp+1,ntortyp-1
1972           nterm_kcc(j,i)=2
1973           nterm_kcc_Tb(j,i)=3
1974           do k=1,nterm_kcc_Tb(j,i)
1975             do l=1,nterm_kcc_Tb(j,i)
1976               v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
1977                               +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1978               v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
1979                               +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1980             enddo
1981           enddo
1982           do k=1,nterm_kcc_Tb(j,i)
1983             do l=1,nterm_kcc_Tb(j,i)
1984 #ifdef CORRCD
1985               v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
1986                               -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1987               v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
1988                               +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1989 #else
1990               v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
1991                               -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1992               v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
1993                               +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1994 #endif
1995             enddo
1996           enddo
1997 !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)
1998         enddo
1999       enddo
2000 #else
2001       write (iout,*) "TOR_MODE>1 only with NEWCORR"
2002       stop
2003 #endif
2004       ENDIF ! TOR_MODE
2005       if (tor_mode.gt.0 .and. lprint) then
2006 !c Print valence-torsional parameters
2007         write (iout,'(a)') &
2008          "Parameters of the valence-torsional potentials"
2009         do i=-ntortyp+1,ntortyp-1
2010         do j=-ntortyp+1,ntortyp-1
2011         write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
2012         write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
2013         do k=1,nterm_kcc(j,i)
2014           do l=1,nterm_kcc_Tb(j,i)
2015             do ll=1,nterm_kcc_Tb(j,i)
2016                write (iout,'(3i5,2f15.4)')&
2017                  k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2018             enddo
2019           enddo
2020         enddo
2021         enddo
2022         enddo
2023       endif
2024
2025 #endif
2026 !elwrite(iout,*) "parmread kontrol sc-bb"
2027 ! Read of Side-chain backbone correlation parameters
2028 ! Modified 11 May 2012 by Adasko
2029 !CC
2030 !
2031      read (isccor,*) nsccortyp
2032
2033      maxinter=3
2034 !c maxinter is maximum interaction sites
2035 !write(iout,*)"maxterm_sccor",maxterm_sccor
2036 !el from module energy-------------
2037       allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2038       allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
2039       allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
2040       allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp))   !(maxterm_sccor,20,20)
2041 !-----------------------------------
2042       allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
2043 !-----------------------------------
2044       allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2045       allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2046                -nsccortyp:nsccortyp))
2047       allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2048                -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2049       allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
2050                -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2051 !-----------------------------------
2052       do i=-nsccortyp,nsccortyp
2053         do j=-nsccortyp,nsccortyp
2054           nterm_sccor(j,i)=0
2055         enddo
2056       enddo
2057 !-----------------------------------
2058
2059       read (isccor,*) (isccortyp(i),i=1,ntyp)
2060       do i=-ntyp,-1
2061         isccortyp(i)=-isccortyp(-i)
2062       enddo
2063       iscprol=isccortyp(20)
2064 !      write (iout,*) 'ntortyp',ntortyp
2065 !      maxinter=3
2066 !c maxinter is maximum interaction sites
2067       do l=1,maxinter
2068       do i=1,nsccortyp
2069         do j=1,nsccortyp
2070           read (isccor,*) &
2071       nterm_sccor(i,j),nlor_sccor(i,j)
2072           v0ijsccor=0.0d0
2073           v0ijsccor1=0.0d0
2074           v0ijsccor2=0.0d0
2075           v0ijsccor3=0.0d0
2076           si=-1.0d0
2077           nterm_sccor(-i,j)=nterm_sccor(i,j)
2078           nterm_sccor(-i,-j)=nterm_sccor(i,j)
2079           nterm_sccor(i,-j)=nterm_sccor(i,j)
2080           do k=1,nterm_sccor(i,j)
2081             read (isccor,*) kk,v1sccor(k,l,i,j),&
2082             v2sccor(k,l,i,j)
2083             if (j.eq.iscprol) then
2084              if (i.eq.isccortyp(10)) then
2085              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2086              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2087              else
2088              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2089                               +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2090              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2091                               +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2092              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2093              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2094              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2095              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2096              endif
2097             else
2098              if (i.eq.isccortyp(10)) then
2099              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2100              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2101              else
2102                if (j.eq.isccortyp(10)) then
2103              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2104              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2105                else
2106              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2107              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2108              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2109              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2110              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2111              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2112                 endif
2113                endif
2114             endif
2115             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2116             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2117             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2118             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2119             si=-si
2120            enddo
2121           do k=1,nlor_sccor(i,j)
2122             read (isccor,*) kk,vlor1sccor(k,i,j),&
2123               vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2124             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2125       (1+vlor3sccor(k,i,j)**2)
2126           enddo
2127           v0sccor(l,i,j)=v0ijsccor
2128           v0sccor(l,-i,j)=v0ijsccor1
2129           v0sccor(l,i,-j)=v0ijsccor2
2130           v0sccor(l,-i,-j)=v0ijsccor3
2131           enddo
2132         enddo
2133       enddo
2134       close (isccor)
2135       if (lprint) then
2136         write (iout,'(/a/)') 'Torsional constants of SCCORR:'
2137         do i=1,nsccortyp
2138           do j=1,nsccortyp
2139             write (iout,*) 'ityp',i,' jtyp',j
2140             write (iout,*) 'Fourier constants'
2141             do k=1,nterm_sccor(i,j)
2142               write (iout,'(2(1pe15.5))') &
2143          (v1sccor(k,l,i,j),v2sccor(k,l,i,j),l=1,maxinter)
2144             enddo
2145             write (iout,*) 'Lorenz constants'
2146             do k=1,nlor_sccor(i,j)
2147               write (iout,'(3(1pe15.5))') &
2148                vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2149             enddo
2150           enddo
2151         enddo
2152       endif
2153
2154 !
2155 ! Read electrostatic-interaction parameters
2156 !
2157       if (lprint) then
2158         write (iout,'(/a)') 'Electrostatic interaction constants:'
2159         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2160                   'IT','JT','APP','BPP','AEL6','AEL3'
2161       endif
2162       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
2163       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
2164       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
2165       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
2166       close (ielep)
2167       do i=1,2
2168         do j=1,2
2169         rri=rpp(i,j)**6
2170         app (i,j)=epp(i,j)*rri*rri
2171         bpp (i,j)=-2.0D0*epp(i,j)*rri
2172         ael6(i,j)=elpp6(i,j)*4.2D0**6
2173         ael3(i,j)=elpp3(i,j)*4.2D0**3
2174         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2175                           ael6(i,j),ael3(i,j)
2176         enddo
2177       enddo
2178 !
2179 ! Read side-chain interaction parameters.
2180 !
2181 !el from module energy - COMMON.INTERACT-------
2182       allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2183       allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2184       allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2185       allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2186       allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2187       allocate(epslip(ntyp,ntyp))
2188       do i=1,ntyp
2189         do j=1,ntyp
2190           augm(i,j)=0.0D0
2191         enddo
2192         chip(i)=0.0D0
2193         alp(i)=0.0D0
2194         sigma0(i)=0.0D0
2195         sigii(i)=0.0D0
2196         rr0(i)=0.0D0
2197       enddo
2198 !--------------------------------
2199
2200       read (isidep,*) ipot,expon
2201 !el      if (ipot.lt.1 .or. ipot.gt.5) then
2202 !        write (iout,'(2a)') 'Error while reading SC interaction',&
2203 !                     'potential file - unknown potential type.'
2204 !        stop
2205 !wl      endif
2206       expon2=expon/2
2207       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2208        ', exponents are ',expon,2*expon
2209 !      goto (10,20,30,30,40) ipot
2210       select case(ipot)
2211 !----------------------- LJ potential ---------------------------------
2212        case (1)
2213 !   10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
2214         read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
2215         if (lprint) then
2216           write (iout,'(/a/)') 'Parameters of the LJ potential:'
2217           write (iout,'(a/)') 'The epsilon array:'
2218           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2219           write (iout,'(/a)') 'One-body parameters:'
2220           write (iout,'(a,4x,a)') 'residue','sigma'
2221           write (iout,'(a3,6x,f10.5)') (restyp(i,molnum(i)),sigma0(i),i=1,ntyp)
2222         endif
2223 !      goto 50
2224 !----------------------- LJK potential --------------------------------
2225        case (2)
2226 !   20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2227         read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2228           (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2229         if (lprint) then
2230           write (iout,'(/a/)') 'Parameters of the LJK potential:'
2231           write (iout,'(a/)') 'The epsilon array:'
2232           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2233           write (iout,'(/a)') 'One-body parameters:'
2234           write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
2235           write (iout,'(a3,6x,2f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),&
2236                 i=1,ntyp)
2237         endif
2238 !      goto 50
2239 !---------------------- GB or BP potential -----------------------------
2240        case (3:4)
2241 !   30 do i=1,ntyp
2242         if (scelemode.eq.0) then
2243         do i=1,ntyp
2244          read (isidep,*)(eps(i,j),j=i,ntyp)
2245         enddo
2246         read (isidep,*)(sigma0(i),i=1,ntyp)
2247         read (isidep,*)(sigii(i),i=1,ntyp)
2248         read (isidep,*)(chip(i),i=1,ntyp)
2249         read (isidep,*)(alp(i),i=1,ntyp)
2250         do i=1,ntyp
2251          read (isidep,*)(epslip(i,j),j=i,ntyp)
2252         enddo
2253 ! For the GB potential convert sigma'**2 into chi'
2254         if (ipot.eq.4) then
2255           do i=1,ntyp
2256             chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2257           enddo
2258         endif
2259         if (lprint) then
2260           write (iout,'(/a/)') 'Parameters of the BP potential:'
2261           write (iout,'(a/)') 'The epsilon array:'
2262           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2263           write (iout,'(/a)') 'One-body parameters:'
2264           write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',&
2265                '    chip  ','    alph  '
2266           write (iout,'(a3,6x,4f10.5)') (restyp(i,molnum(i)),sigma0(i),sigii(i),&
2267                            chip(i),alp(i),i=1,ntyp)
2268         endif
2269         else
2270       allocate(icharge(ntyp1))
2271 !      print *,ntyp,icharge(i)
2272       icharge(:)=0
2273       read (isidep,*) (icharge(i),i=1,ntyp)
2274       print *,ntyp,icharge(i)
2275 !      if(.not.allocated(eps)) allocate(eps(-ntyp
2276       write (2,*) "icharge",(icharge(i),i=1,ntyp)
2277        allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2278        allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2279        allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2280        allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2281        allocate(epsintab(ntyp,ntyp))
2282        allocate(dtail(2,ntyp,ntyp))
2283       print *,"control line 1"
2284        allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2285        allocate(wqdip(2,ntyp,ntyp))
2286        allocate(wstate(4,ntyp,ntyp))
2287        allocate(dhead(2,2,ntyp,ntyp))
2288        allocate(nstate(ntyp,ntyp))
2289        allocate(debaykap(ntyp,ntyp))
2290       print *,"control line 2"
2291       if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2292       if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2293
2294       do i=1,ntyp
2295        do j=1,i
2296 !        write (*,*) "Im in ALAB", i, " ", j
2297         read(isidep,*) &
2298        eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), &
2299        (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &
2300        chis(i,j),chis(j,i), &
2301        nstate(i,j),(wstate(k,i,j),k=1,4), &
2302        dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),&
2303        dtail(1,i,j),dtail(2,i,j), &
2304        epshead(i,j),sig0head(i,j), &
2305        rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), &
2306        alphapol(i,j),alphapol(j,i), &
2307        (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j)
2308 !       print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
2309        END DO
2310       END DO
2311       DO i = 1, ntyp
2312        DO j = i+1, ntyp
2313         eps(i,j) = eps(j,i)
2314         sigma(i,j) = sigma(j,i)
2315         sigmap1(i,j)=sigmap1(j,i)
2316         sigmap2(i,j)=sigmap2(j,i)
2317         sigiso1(i,j)=sigiso1(j,i)
2318         sigiso2(i,j)=sigiso2(j,i)
2319 !        print *,"ATU",sigma(j,i),sigma(i,j),i,j
2320         nstate(i,j) = nstate(j,i)
2321         dtail(1,i,j) = dtail(1,j,i)
2322         dtail(2,i,j) = dtail(2,j,i)
2323         DO k = 1, 4
2324          alphasur(k,i,j) = alphasur(k,j,i)
2325          wstate(k,i,j) = wstate(k,j,i)
2326          alphiso(k,i,j) = alphiso(k,j,i)
2327         END DO
2328
2329         dhead(2,1,i,j) = dhead(1,1,j,i)
2330         dhead(2,2,i,j) = dhead(1,2,j,i)
2331         dhead(1,1,i,j) = dhead(2,1,j,i)
2332         dhead(1,2,i,j) = dhead(2,2,j,i)
2333
2334         epshead(i,j) = epshead(j,i)
2335         sig0head(i,j) = sig0head(j,i)
2336
2337         DO k = 1, 2
2338          wqdip(k,i,j) = wqdip(k,j,i)
2339         END DO
2340
2341         wquad(i,j) = wquad(j,i)
2342         epsintab(i,j) = epsintab(j,i)
2343         debaykap(i,j)=debaykap(j,i)
2344 !        if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2345        END DO
2346       END DO
2347       endif
2348
2349 !      goto 50
2350 !--------------------- GBV potential -----------------------------------
2351        case (5)
2352 !   40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2353         read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2354           (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2355         (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2356         if (lprint) then
2357           write (iout,'(/a/)') 'Parameters of the GBV potential:'
2358           write (iout,'(a/)') 'The epsilon array:'
2359           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2360           write (iout,'(/a)') 'One-body parameters:'
2361           write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',&
2362             's||/s_|_^2','    chip  ','    alph  '
2363           write (iout,'(a3,6x,5f10.5)') (restyp(i,molnum(i)),sigma0(i),rr0(i),&
2364                  sigii(i),chip(i),alp(i),i=1,ntyp)
2365         endif
2366        case default
2367         write (iout,'(2a)') 'Error while reading SC interaction',&
2368                      'potential file - unknown potential type.'
2369         stop
2370 !   50 continue
2371       end select
2372 !      continue
2373       close (isidep)
2374 !-----------------------------------------------------------------------
2375 ! Calculate the "working" parameters of SC interactions.
2376
2377 !el from module energy - COMMON.INTERACT-------
2378 !      allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2379             if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1))
2380       allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1)) !(ntyp,ntyp)
2381       allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2382       if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2383       allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2384       allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2385         dcavtub(ntyp1))
2386       allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2387         tubetranene(ntyp1))
2388       do i=1,ntyp1
2389         do j=1,ntyp1
2390           aa_aq(i,j)=0.0D0
2391           bb_aq(i,j)=0.0D0
2392           aa_lip(i,j)=0.0d0
2393           bb_lip(i,j)=0.0d0
2394                if (scelemode.eq.0) then
2395           chi(i,j)=0.0D0
2396           sigma(i,j)=0.0D0
2397           r0(i,j)=0.0D0
2398             endif
2399         enddo
2400       enddo
2401 !--------------------------------
2402
2403       do i=2,ntyp
2404         do j=1,i-1
2405           eps(i,j)=eps(j,i)
2406           epslip(i,j)=epslip(j,i)
2407         enddo
2408       enddo
2409       if (scelemode.eq.0) then
2410       do i=1,ntyp
2411         do j=i,ntyp
2412           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2413           sigma(j,i)=sigma(i,j)
2414           rs0(i,j)=dwa16*sigma(i,j)
2415           rs0(j,i)=rs0(i,j)
2416         enddo
2417       enddo
2418       endif
2419       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2420        'Working parameters of the SC interactions:',&
2421        '     a    ','     b    ','   augm   ','  sigma ','   r0   ',&
2422        '  chi1   ','   chi2   '
2423       do i=1,ntyp
2424         do j=i,ntyp
2425           epsij=eps(i,j)
2426           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2427             rrij=sigma(i,j)
2428             print *,"rrij",rrij
2429           else
2430             rrij=rr0(i)+rr0(j)
2431           endif
2432           r0(i,j)=rrij
2433           r0(j,i)=rrij
2434           rrij=rrij**expon
2435           epsij=eps(i,j)
2436           sigeps=dsign(1.0D0,epsij)
2437           epsij=dabs(epsij)
2438           aa_aq(i,j)=epsij*rrij*rrij
2439           bb_aq(i,j)=-sigeps*epsij*rrij
2440           aa_aq(j,i)=aa_aq(i,j)
2441           bb_aq(j,i)=bb_aq(i,j)
2442           epsijlip=epslip(i,j)
2443           sigeps=dsign(1.0D0,epsijlip)
2444           epsijlip=dabs(epsijlip)
2445           aa_lip(i,j)=epsijlip*rrij*rrij
2446           bb_lip(i,j)=-sigeps*epsijlip*rrij
2447           aa_lip(j,i)=aa_lip(i,j)
2448           bb_lip(j,i)=bb_lip(i,j)
2449           if ((ipot.gt.2).and. (scelemode.eq.0))then
2450             sigt1sq=sigma0(i)**2
2451             sigt2sq=sigma0(j)**2
2452             sigii1=sigii(i)
2453             sigii2=sigii(j)
2454             ratsig1=sigt2sq/sigt1sq
2455             ratsig2=1.0D0/ratsig1
2456             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2457             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2458             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2459           else
2460             rsum_max=sigma(i,j)
2461           endif
2462 !         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2463             sigmaii(i,j)=rsum_max
2464             sigmaii(j,i)=rsum_max
2465 !         else
2466 !           sigmaii(i,j)=r0(i,j)
2467 !           sigmaii(j,i)=r0(i,j)
2468 !         endif
2469 !d        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2470           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2471             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2472             augm(i,j)=epsij*r_augm**(2*expon)
2473 !           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2474             augm(j,i)=augm(i,j)
2475           else
2476             augm(i,j)=0.0D0
2477             augm(j,i)=0.0D0
2478           endif
2479           if (lprint) then
2480             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')  &
2481             restyp(i,molnum(i)),restyp(j,molnum(i)),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2482             sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2483           endif
2484         enddo
2485       enddo
2486       allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2487       allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2488       allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2489       allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2490       allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2491       allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2492       allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2493       allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2494       allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2495       allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2496       allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2497       allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2498       allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2499       allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2500       allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2501       allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2502
2503 !      augm(:,:)=0.0D0
2504 !      chip(:)=0.0D0
2505 !      alp(:)=0.0D0
2506 !      sigma0(:)=0.0D0
2507 !      sigii(:)=0.0D0
2508 !      rr0(:)=0.0D0
2509
2510       read (isidep_nucl,*) ipot_nucl
2511 !      print *,"TU?!",ipot_nucl
2512       if (ipot_nucl.eq.1) then
2513         do i=1,ntyp_molec(2)
2514           do j=i,ntyp_molec(2)
2515             read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2516             elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2517           enddo
2518         enddo
2519       else
2520         do i=1,ntyp_molec(2)
2521           do j=i,ntyp_molec(2)
2522             read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2523                chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2524                elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2525           enddo
2526         enddo
2527       endif
2528 !      rpp(1,1)=2**(1.0/6.0)*5.16158
2529       do i=1,ntyp_molec(2)
2530         do j=i,ntyp_molec(2)
2531           rrij=sigma_nucl(i,j)
2532           r0_nucl(i,j)=rrij
2533           r0_nucl(j,i)=rrij
2534           rrij=rrij**expon
2535           epsij=4*eps_nucl(i,j)
2536           sigeps=dsign(1.0D0,epsij)
2537           epsij=dabs(epsij)
2538           aa_nucl(i,j)=epsij*rrij*rrij
2539           bb_nucl(i,j)=-sigeps*epsij*rrij
2540           ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2541           ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2542           ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2543           ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2544           sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2545          2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2546         enddo
2547         do j=1,i-1
2548           aa_nucl(i,j)=aa_nucl(j,i)
2549           bb_nucl(i,j)=bb_nucl(j,i)
2550           ael3_nucl(i,j)=ael3_nucl(j,i)
2551           ael6_nucl(i,j)=ael6_nucl(j,i)
2552           ael63_nucl(i,j)=ael63_nucl(j,i)
2553           ael32_nucl(i,j)=ael32_nucl(j,i)
2554           elpp3_nucl(i,j)=elpp3_nucl(j,i)
2555           elpp6_nucl(i,j)=elpp6_nucl(j,i)
2556           elpp63_nucl(i,j)=elpp63_nucl(j,i)
2557           elpp32_nucl(i,j)=elpp32_nucl(j,i)
2558           eps_nucl(i,j)=eps_nucl(j,i)
2559           sigma_nucl(i,j)=sigma_nucl(j,i)
2560           sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2561         enddo
2562       enddo
2563
2564       write(iout,*) "tube param"
2565       read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2566       ccavtubpep,dcavtubpep,tubetranenepep
2567       sigmapeptube=sigmapeptube**6
2568       sigeps=dsign(1.0D0,epspeptube)
2569       epspeptube=dabs(epspeptube)
2570       pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2571       pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2572       write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2573       do i=1,ntyp
2574        read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2575       ccavtub(i),dcavtub(i),tubetranene(i)
2576        sigmasctube=sigmasctube**6
2577        sigeps=dsign(1.0D0,epssctube)
2578        epssctube=dabs(epssctube)
2579        sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2580        sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2581       write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2582       enddo
2583 !-----------------READING SC BASE POTENTIALS-----------------------------
2584       allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2585       allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2586       allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2587       allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2588       allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2589       allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2590       allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2591       allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2592       allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2593       allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2594       allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2595       allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2596       allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2597       allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2598       allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2599       allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2600
2601
2602       do i=1,ntyp_molec(1)
2603        do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2604         write (*,*) "Im in ", i, " ", j
2605         read(isidep_scbase,*) &
2606         eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2607         chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2608          write(*,*) "eps",eps_scbase(i,j)
2609         read(isidep_scbase,*) &
2610        (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2611        chis_scbase(i,j,1),chis_scbase(i,j,2)
2612         read(isidep_scbase,*) &
2613        dhead_scbasei(i,j), &
2614        dhead_scbasej(i,j), &
2615        rborn_scbasei(i,j),rborn_scbasej(i,j)
2616         read(isidep_scbase,*) &
2617        (wdipdip_scbase(k,i,j),k=1,3), &
2618        (wqdip_scbase(k,i,j),k=1,2)
2619         read(isidep_scbase,*) &
2620        alphapol_scbase(i,j), &
2621        epsintab_scbase(i,j)
2622        END DO
2623       END DO
2624       allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2625       allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2626
2627       do i=1,ntyp_molec(1)
2628        do j=1,ntyp_molec(2)-1
2629           epsij=eps_scbase(i,j)
2630           rrij=sigma_scbase(i,j)
2631 !          r0(i,j)=rrij
2632 !          r0(j,i)=rrij
2633           rrij=rrij**expon
2634 !          epsij=eps(i,j)
2635           sigeps=dsign(1.0D0,epsij)
2636           epsij=dabs(epsij)
2637           aa_scbase(i,j)=epsij*rrij*rrij
2638           bb_scbase(i,j)=-sigeps*epsij*rrij
2639         enddo
2640        enddo
2641 !-----------------READING PEP BASE POTENTIALS-------------------
2642       allocate(eps_pepbase(ntyp_molec(2)))
2643       allocate(sigma_pepbase(ntyp_molec(2)))
2644       allocate(chi_pepbase(ntyp_molec(2),2))
2645       allocate(chipp_pepbase(ntyp_molec(2),2))
2646       allocate(alphasur_pepbase(4,ntyp_molec(2)))
2647       allocate(sigmap1_pepbase(ntyp_molec(2)))
2648       allocate(sigmap2_pepbase(ntyp_molec(2)))
2649       allocate(chis_pepbase(ntyp_molec(2),2))
2650       allocate(wdipdip_pepbase(3,ntyp_molec(2)))
2651
2652
2653        do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2654         write (*,*) "Im in ", i, " ", j
2655         read(isidep_pepbase,*) &
2656         eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
2657         chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
2658          write(*,*) "eps",eps_pepbase(j)
2659         read(isidep_pepbase,*) &
2660        (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
2661        chis_pepbase(j,1),chis_pepbase(j,2)
2662         read(isidep_pepbase,*) &
2663        (wdipdip_pepbase(k,j),k=1,3)
2664        END DO
2665       allocate(aa_pepbase(ntyp_molec(2)))
2666       allocate(bb_pepbase(ntyp_molec(2)))
2667
2668        do j=1,ntyp_molec(2)-1
2669           epsij=eps_pepbase(j)
2670           rrij=sigma_pepbase(j)
2671 !          r0(i,j)=rrij
2672 !          r0(j,i)=rrij
2673           rrij=rrij**expon
2674 !          epsij=eps(i,j)
2675           sigeps=dsign(1.0D0,epsij)
2676           epsij=dabs(epsij)
2677           aa_pepbase(j)=epsij*rrij*rrij
2678           bb_pepbase(j)=-sigeps*epsij*rrij
2679         enddo
2680 !--------------READING SC PHOSPHATE-------------------------------------
2681 !--------------READING SC PHOSPHATE-------------------------------------
2682       allocate(eps_scpho(ntyp_molec(1)))
2683       allocate(sigma_scpho(ntyp_molec(1)))
2684       allocate(chi_scpho(ntyp_molec(1),2))
2685       allocate(chipp_scpho(ntyp_molec(1),2))
2686       allocate(alphasur_scpho(4,ntyp_molec(1)))
2687       allocate(sigmap1_scpho(ntyp_molec(1)))
2688       allocate(sigmap2_scpho(ntyp_molec(1)))
2689       allocate(chis_scpho(ntyp_molec(1),2))
2690       allocate(wqq_scpho(ntyp_molec(1)))
2691       allocate(wqdip_scpho(2,ntyp_molec(1)))
2692       allocate(alphapol_scpho(ntyp_molec(1)))
2693       allocate(epsintab_scpho(ntyp_molec(1)))
2694       allocate(dhead_scphoi(ntyp_molec(1)))
2695       allocate(rborn_scphoi(ntyp_molec(1)))
2696       allocate(rborn_scphoj(ntyp_molec(1)))
2697       allocate(alphi_scpho(ntyp_molec(1)))
2698
2699
2700 !      j=1
2701        do j=1,ntyp_molec(1) ! without U then we will take T for U
2702         write (*,*) "Im in scpho ", i, " ", j
2703         read(isidep_scpho,*) &
2704         eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
2705         chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
2706          write(*,*) "eps",eps_scpho(j)
2707         read(isidep_scpho,*) &
2708        (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
2709        chis_scpho(j,1),chis_scpho(j,2)
2710         read(isidep_scpho,*) &
2711        (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
2712         read(isidep_scpho,*) &
2713          epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
2714          alphi_scpho(j)
2715
2716        END DO
2717       allocate(aa_scpho(ntyp_molec(1)))
2718       allocate(bb_scpho(ntyp_molec(1)))
2719        do j=1,ntyp_molec(1)
2720           epsij=eps_scpho(j)
2721           rrij=sigma_scpho(j)
2722 !          r0(i,j)=rrij
2723 !          r0(j,i)=rrij
2724           rrij=rrij**expon
2725 !          epsij=eps(i,j)
2726           sigeps=dsign(1.0D0,epsij)
2727           epsij=dabs(epsij)
2728           aa_scpho(j)=epsij*rrij*rrij
2729           bb_scpho(j)=-sigeps*epsij*rrij
2730         enddo
2731
2732
2733         read(isidep_peppho,*) &
2734         eps_peppho,sigma_peppho
2735         read(isidep_peppho,*) &
2736        (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
2737         read(isidep_peppho,*) &
2738        (wqdip_peppho(k),k=1,2)
2739
2740           epsij=eps_peppho
2741           rrij=sigma_peppho
2742 !          r0(i,j)=rrij
2743 !          r0(j,i)=rrij
2744           rrij=rrij**expon
2745 !          epsij=eps(i,j)
2746           sigeps=dsign(1.0D0,epsij)
2747           epsij=dabs(epsij)
2748           aa_peppho=epsij*rrij*rrij
2749           bb_peppho=-sigeps*epsij*rrij
2750
2751
2752       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2753       do i=1,ntyp
2754         do j=1,2
2755           bad(i,j)=0.0D0
2756         enddo
2757       enddo
2758 ! Ions by Aga
2759
2760        allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
2761        allocate(alphapolcat2(ntyp,ntyp))
2762        allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
2763        allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
2764        allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp))
2765        allocate(epsintabcat(ntyp,ntyp))
2766        allocate(dtailcat(2,ntyp,ntyp))
2767        allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
2768        allocate(wqdipcat(2,ntyp,ntyp))
2769        allocate(wstatecat(4,ntyp,ntyp))
2770        allocate(dheadcat(2,2,ntyp,ntyp))
2771        allocate(nstatecat(ntyp,ntyp))
2772        allocate(debaykapcat(ntyp,ntyp))
2773
2774       if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1))
2775       if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
2776 !      if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
2777       if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp)
2778       if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
2779
2780
2781       allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5)))
2782 ! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
2783        if (oldion.eq.0) then
2784             if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
2785             allocate(icharge(1:ntyp1))
2786             read(iion,*) (icharge(i),i=1,ntyp)
2787             else
2788              read(iion,*) ijunk
2789             endif
2790
2791             do i=-ntyp_molec(5),ntyp_molec(5)
2792              read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
2793              print *,msc(i,5),restok(i,5)
2794             enddo
2795             ip(5)=0.2
2796 !DIR$ NOUNROLL
2797       do j=1,ntyp_molec(5)-1
2798        do i=1,ntyp
2799 !       do j=1,ntyp_molec(5)
2800 !        write (*,*) "Im in ALAB", i, " ", j
2801         read(iion,*) &
2802        epscat(i,j),sigmacat(i,j), &
2803 !       chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
2804        chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), &
2805
2806        (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
2807 !       chiscat(i,j),chiscat(j,i), &
2808        chis1cat(i,j),chis2cat(i,j), &
2809
2810        nstatecat(i,j),(wstatecat(k,i,j),k=1,4), &                           !5 w tej lini - 1 integer pierwszy
2811        dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
2812        dtailcat(1,i,j),dtailcat(2,i,j), &
2813        epsheadcat(i,j),sig0headcat(i,j), &
2814 !wdipcat = w1 , w2
2815 !       rborncat(i,j),rborncat(j,i),&
2816        rborn1cat(i,j),rborn2cat(i,j),&
2817        (wqdipcat(k,i,j),k=1,2), &
2818        alphapolcat(i,j),alphapolcat2(j,i), &
2819        (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
2820 !       print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
2821        END DO
2822       END DO
2823       allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
2824       do i=1,ntyp
2825         do j=1,ntyp_molec(5)-1 !without zinc
2826           epsij=epscat(i,j)
2827           rrij=sigmacat(i,j)
2828           rrij=rrij**expon
2829           sigeps=dsign(1.0D0,epsij)
2830           epsij=dabs(epsij)
2831           aa_aq_cat(i,j)=epsij*rrij*rrij
2832           bb_aq_cat(i,j)=-sigeps*epsij*rrij
2833          enddo
2834        enddo
2835        do i=1,ntyp
2836        do j=1,ntyp_molec(5)
2837       if (i.eq.10) then
2838       write (iout,*) 'i= ', i, ' j= ', j
2839       write (iout,*) 'epsi0= ', epscat(i,j)
2840       write (iout,*) 'sigma0= ', sigmacat(i,j)
2841       write (iout,*) 'chi1= ', chi1cat(i,j)
2842       write (iout,*) 'chi1= ', chi2cat(i,j)
2843       write (iout,*) 'chip1= ', chipp1cat(1,j)
2844       write (iout,*) 'chip2= ', chipp2cat(1,j)
2845       write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
2846       write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
2847       write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
2848       write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
2849       write (iout,*) 'sig1= ', sigmap1cat(1,j)
2850       write (iout,*) 'sig2= ', sigmap2cat(1,j)
2851       write (iout,*) 'chis1= ', chis1cat(1,j)
2852       write (iout,*) 'chis1= ', chis2cat(1,j)
2853       write (iout,*) 'nstatecat(i,j)= ', nstatecat(1,j)
2854       write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,1,j)
2855       write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
2856       write (iout,*) 'dhead2= ', dheadcat(1,2,1,j)
2857       write (iout,*) 'a1= ', rborn1cat(i,j)
2858       write (iout,*) 'a2= ', rborn2cat(i,j)
2859       write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
2860       write (iout,*) 'alphapol1= ',  alphapolcat(1,j)
2861       write (iout,*) 'alphapol2= ',  alphapolcat(j,1)
2862       write (iout,*) 'w1= ', wqdipcat(1,i,j)
2863       write (iout,*) 'w2= ', wqdipcat(2,i,j)
2864       write (iout,*) 'debaykapcat(i,j)= ',  debaykapcat(1,j)
2865       endif
2866
2867       If ((i.eq.1).and.(j.eq.27)) then
2868       write (iout,*) 'SEP'
2869       Write (iout,*) 'w1= ', wqdipcat(1,1,27)
2870       Write (iout,*) 'w2= ', wqdipcat(2,1,27)
2871       endif
2872
2873        enddo
2874        enddo
2875
2876       endif
2877 ! here two denotes the Zn2+ and Cu2+
2878       write(iout,*) "before TRANPARM"
2879       allocate(aomicattr(0:3,2))
2880       allocate(athetacattran(0:6,5,2))
2881       allocate(agamacattran(3,5,2))
2882       allocate(acatshiftdsc(5,2))
2883       allocate(bcatshiftdsc(5,2))
2884       allocate(demorsecat(5,2))
2885       allocate(alphamorsecat(5,2))
2886       allocate(x0catleft(5,2))
2887       allocate(x0catright(5,2))
2888       allocate(x0cattrans(5,2))
2889       allocate(ntrantyp(2))
2890       do i=1,1 ! currently only Zn2+
2891
2892       read(iiontran,*) ntrantyp(i)
2893 !ntrantyp=4
2894 !| ao0 ao1 ao2 ao3
2895 !ASP| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
2896 !CYS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0left x0right x0transi
2897 !GLU| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -0.5
2898 !HIS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
2899       read(iiontran,*) (aomicattr(j,i),j=0,3)
2900       do j=1,ntrantyp(i)
2901        read (iiontran,*) (agamacattran(k,j,i),k=1,3),&
2902        (athetacattran(k,j,i),k=0,6),acatshiftdsc(j,i),bcatshiftdsc(j,i),&
2903        demorsecat(j,i),alphamorsecat(j,i),x0catleft(j,i),x0catright(j,i),&
2904        x0cattrans(j,i)
2905       enddo
2906       enddo
2907 #ifdef CLUSTER
2908 !
2909 ! Define the SC-p interaction constants
2910 !
2911       do i=1,20
2912         do j=1,2
2913           eps_scp(i,j)=-1.5d0
2914           rscp(i,j)=4.0d0
2915         enddo
2916       enddo
2917 #endif
2918       allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2919       read (itorp_nucl,*) ntortyp_nucl
2920 !      print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2921 !el from energy module---------
2922       allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2923       allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2924
2925       allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2926       allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2927       allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2928       allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2929
2930       allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2931       allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2932 !el---------------------------
2933       nterm_nucl(:,:)=0
2934       nlor_nucl(:,:)=0
2935 !el--------------------
2936       read (itorp_nucl,*) &
2937         (itortyp_nucl(i),i=1,ntyp_molec(2))
2938 !        print *,itortyp_nucl(:)
2939 !c      write (iout,*) 'ntortyp',ntortyp
2940       do i=1,ntortyp_nucl
2941         do j=1,ntortyp_nucl
2942           read (itorp_nucl,*) nterm_nucl(i,j),nlor_nucl(i,j)
2943 !           print *,nterm_nucl(i,j),nlor_nucl(i,j)
2944           v0ij=0.0d0
2945           si=-1.0d0
2946           do k=1,nterm_nucl(i,j)
2947             read (itorp_nucl,*) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2948             v0ij=v0ij+si*v1_nucl(k,i,j)
2949             si=-si
2950           enddo
2951           do k=1,nlor_nucl(i,j)
2952             read (itorp_nucl,*) kk,vlor1_nucl(k,i,j),&
2953               vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2954             v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2955           enddo
2956           v0_nucl(i,j)=v0ij
2957         enddo
2958       enddo
2959
2960
2961 !elwrite(iout,*) "parmread kontrol before oldscp"
2962 !
2963 ! Define the SC-p interaction constants
2964 !
2965 #ifdef OLDSCP
2966       do i=1,20
2967 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2968 ! helix formation)
2969 !       aad(i,1)=0.3D0*4.0D0**12
2970 ! Following line for constants currently implemented
2971 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2972         aad(i,1)=1.5D0*4.0D0**12
2973 !       aad(i,1)=0.17D0*5.6D0**12
2974         aad(i,2)=aad(i,1)
2975 ! "Soft" SC-p repulsion
2976         bad(i,1)=0.0D0
2977 ! Following line for constants currently implemented
2978 !       aad(i,1)=0.3D0*4.0D0**6
2979 ! "Hard" SC-p repulsion
2980         bad(i,1)=3.0D0*4.0D0**6
2981 !       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2982         bad(i,2)=bad(i,1)
2983 !       aad(i,1)=0.0D0
2984 !       aad(i,2)=0.0D0
2985 !       bad(i,1)=1228.8D0
2986 !       bad(i,2)=1228.8D0
2987       enddo
2988 #else
2989 !
2990 ! 8/9/01 Read the SC-p interaction constants from file
2991 !
2992       do i=1,ntyp
2993         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
2994       enddo
2995       do i=1,ntyp
2996         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2997         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2998         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2999         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
3000       enddo
3001
3002       if (lprint) then
3003         write (iout,*) "Parameters of SC-p interactions:"
3004         do i=1,20
3005           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
3006            eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
3007         enddo
3008       endif
3009 #endif
3010       allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
3011
3012       do i=1,ntyp_molec(2)
3013         read (iscpp_nucl,*) eps_scp_nucl(i),rscp_nucl(i)
3014       enddo
3015       do i=1,ntyp_molec(2)
3016         aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
3017         bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
3018       enddo
3019       r0pp=1.12246204830937298142*5.16158
3020       epspp=4.95713/4
3021       AEES=108.661
3022       BEES=0.433246
3023
3024 !
3025 ! Define the constants of the disulfide bridge
3026 !
3027 !      ebr=-5.50D0
3028 !
3029 ! Old arbitrary potential - commented out.
3030 !
3031 !      dbr= 4.20D0
3032 !      fbr= 3.30D0
3033 !
3034 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
3035 ! energy surface of diethyl disulfide.
3036 ! A. Liwo and U. Kozlowska, 11/24/03
3037 !
3038 !      D0CM = 3.78d0
3039 !      AKCM = 15.1d0
3040 !      AKTH = 11.0d0
3041 !      AKCT = 12.0d0
3042 !      V1SS =-1.08d0
3043 !      V2SS = 7.61d0
3044 !      V3SS = 13.7d0
3045 #ifndef CLUSTER
3046
3047       if (dyn_ss) then
3048        ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
3049         Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
3050         akcm=akcm/wsc*ssscale
3051         akth=akth/wsc*ssscale
3052         akct=akct/wsc*ssscale
3053         v1ss=v1ss/wsc*ssscale
3054         v2ss=v2ss/wsc*ssscale
3055         v3ss=v3ss/wsc*ssscale
3056       else
3057         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
3058       endif
3059 #endif
3060       if (lprint) then
3061       write (iout,'(/a)') "Disulfide bridge parameters:"
3062       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
3063       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
3064       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
3065       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
3066        ' v3ss:',v3ss
3067       endif
3068       return
3069       end subroutine parmread
3070 #ifndef CLUSTER
3071 !-----------------------------------------------------------------------------
3072 ! mygetenv.F
3073 !-----------------------------------------------------------------------------
3074       subroutine mygetenv(string,var)
3075 !
3076 ! Version 1.0
3077 !
3078 ! This subroutine passes the environmental variables to FORTRAN program.
3079 ! If the flags -DMYGETENV and -DMPI are not for compilation, it calls the
3080 ! standard FORTRAN GETENV subroutine. If both flags are set, the subroutine
3081 ! reads the environmental variables from $HOME/.env
3082 !
3083 ! Usage: As for the standard FORTRAN GETENV subroutine.
3084 !
3085 ! Purpose: some versions/installations of MPI do not transfer the environmental
3086 ! variables to slave processors, if these variables are set in the shell script
3087 ! from which mpirun is called.
3088 !
3089 ! A.Liwo, 7/29/01
3090 !
3091 #ifdef MPI
3092       use MPI_data
3093       include "mpif.h"
3094 #endif
3095 !      implicit none
3096       character*(*) :: string,var
3097 #if defined(MYGETENV) && defined(MPI)
3098 !      include "DIMENSIONS.ZSCOPT"
3099 !      include "mpif.h"
3100 !      include "COMMON.MPI"
3101 !el      character*360 ucase
3102 !el      external ucase
3103       character(len=360) :: string1(360),karta
3104       character(len=240) :: home
3105       integer i,n !,ilen
3106 !el      external ilen
3107       call getenv("HOME",home)
3108       open(99,file=home(:ilen(home))//"/.env",status="OLD",err=112)
3109       do while (.true.)
3110         read (99,end=111,err=111,'(a)') karta
3111         do i=1,80
3112           string1(i)=" "
3113         enddo
3114         call split_string(karta,string1,80,n)
3115         if (ucase(string1(1)(:ilen(string1(1)))).eq."SETENV" .and. &
3116          string1(2)(:ilen(string1(2))).eq.string(:ilen(string)) ) then
3117            var=string1(3)
3118            print *,"Processor",me,": ",var(:ilen(var)),&
3119             " assigned to ",string(:ilen(string))
3120            close(99)
3121            return
3122         endif
3123       enddo
3124  111  print *,"Environment variable ",string(:ilen(string))," not set."
3125       close(99)
3126       return
3127  112  print *,"Error opening environment file!"
3128 #else
3129       call getenv(string,var)
3130 #endif
3131       return
3132       end subroutine mygetenv
3133 !-----------------------------------------------------------------------------
3134 ! readrtns.F
3135 !-----------------------------------------------------------------------------
3136       subroutine read_general_data(*)
3137
3138       use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions,&
3139           scelemode,TUBEmode,tor_mode,energy_dec,r_cut_ang,r_cut_mart,&
3140           rlamb_mart
3141
3142       use energy_data, only:distchainmax,tubeR0,tubecenter,dyn_ss,constr_homology
3143       use geometry_data, only:boxxsize,boxysize,boxzsize,bordtubetop,&
3144           bordtubebot,tubebufthick,buftubebot,buftubetop,buflipbot, bufliptop,bordlipbot,bordliptop,     &
3145         lipbufthick,lipthick
3146 !      implicit none
3147 !      include "DIMENSIONS"
3148 !      include "DIMENSIONS.ZSCOPT"
3149 !      include "DIMENSIONS.FREE"
3150 !      include "COMMON.TORSION"
3151 !      include "COMMON.INTERACT"
3152 !      include "COMMON.IOUNITS"
3153 !      include "COMMON.TIME1"
3154 !      include "COMMON.PROT"
3155 !      include "COMMON.PROTFILES"
3156 !      include "COMMON.CHAIN"
3157 !      include "COMMON.NAMES"
3158 !      include "COMMON.FFIELD"
3159 !      include "COMMON.ENEPS"
3160 !      include "COMMON.WEIGHTS"
3161 !      include "COMMON.FREE"
3162 !      include "COMMON.CONTROL"
3163 !      include "COMMON.ENERGIES"
3164
3165       character(len=800) :: controlcard
3166       integer :: i,j,k,ii,n_ene_found
3167       integer :: ind,itype1,itype2,itypf,itypsc,itypp
3168 !el      integer ilen
3169 !el      external ilen
3170 !el      character*16 ucase
3171       character(len=16) :: key
3172 !el      external ucase
3173       call card_concat(controlcard,.true.)
3174       call readi(controlcard,"N_ENE",n_eneW,max_eneW)
3175       if (n_eneW.gt.max_eneW) then
3176         write (iout,*) "Error: parameter out of range: N_ENE",n_eneW,&
3177           max_eneW
3178         return 1
3179       endif
3180       call readi(controlcard,"NPARMSET",nparmset,1)
3181 !elwrite(iout,*)"in read_gen data"
3182       separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
3183       call readi(controlcard,"IPARMPRINT",iparmprint,1)
3184       write (iout,*) "PARMPRINT",iparmprint
3185       if (nparmset.gt.max_parm) then
3186         write (iout,*) "Error: parameter out of range: NPARMSET",&
3187           nparmset, Max_Parm
3188         return 1
3189       endif
3190 !elwrite(iout,*)"in read_gen data"
3191       call readi(controlcard,"MAXIT",maxit,5000)
3192       call reada(controlcard,"FIMIN",fimin,1.0d-3)
3193       call readi(controlcard,"ENSEMBLES",ensembles,0)
3194       hamil_rep=index(controlcard,"HAMIL_REP").gt.0
3195       write (iout,*) "Number of energy parameter sets",nparmset
3196       allocate(isampl(nparmset))
3197       call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
3198       write (iout,*) "MaxSlice",MaxSlice
3199       call readi(controlcard,"NSLICE",nslice,1)
3200 !elwrite(iout,*)"in read_gen data"
3201       call flush(iout)
3202       if (nslice.gt.MaxSlice) then
3203         write (iout,*) "Error: parameter out of range: NSLICE",nslice,&
3204           MaxSlice
3205         return 1
3206       endif
3207       write (iout,*) "Frequency of storing conformations",&
3208        (isampl(i),i=1,nparmset)
3209       write (iout,*) "Maxit",maxit," Fimin",fimin
3210       call readi(controlcard,"NQ",nQ,1)
3211       if (nQ.gt.MaxQ) then
3212         write (iout,*) "Error: parameter out of range: NQ",nq,&
3213           maxq
3214         return 1
3215       endif
3216       indpdb=0
3217       if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
3218       call reada(controlcard,"DELTA",delta,1.0d-2)
3219       call readi(controlcard,"EINICHECK",einicheck,2)
3220       call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
3221       call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
3222       call readi(controlcard,"RESCALE",rescale_modeW,1)
3223       call reada(controlcard,'BOXX',boxxsize,100.0d0)
3224       call reada(controlcard,'BOXY',boxysize,100.0d0)
3225       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
3226       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3227       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3228       if (lipthick.gt.0.0d0) then
3229        bordliptop=(boxzsize+lipthick)/2.0
3230        bordlipbot=bordliptop-lipthick
3231       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3232       write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3233       buflipbot=bordlipbot+lipbufthick
3234       bufliptop=bordliptop-lipbufthick
3235       if ((lipbufthick*2.0d0).gt.lipthick) &
3236        write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3237       endif !lipthick.gt.0
3238       write(iout,*) "bordliptop=",bordliptop
3239       write(iout,*) "bordlipbot=",bordlipbot
3240       write(iout,*) "bufliptop=",bufliptop
3241       write(iout,*) "buflipbot=",buflipbot
3242
3243       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3244       call readi(controlcard,"SCELEMODE",scelemode,0)
3245       call readi(controlcard,"OLDION",oldion,0)
3246       dyn_ss=(index(controlcard,'DYN_SS').gt.0)
3247       print *,"SCELE",scelemode
3248       call readi(controlcard,'TORMODE',tor_mode,0)
3249 !C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3250         write(iout,*) "torsional and valence angle mode",tor_mode
3251
3252       call readi(controlcard,'TUBEMOD',tubemode,0)
3253       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
3254 !c      if (constr_homology) tole=dmax1(tole,1.5d0)
3255       write (iout,*) "with_homology_constr ",with_dihed_constr, &
3256        " CONSTR_HOMOLOGY",constr_homology
3257 !      read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
3258 !      out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
3259 !      out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
3260
3261       if (TUBEmode.gt.0) then
3262        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
3263        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
3264        call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
3265        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
3266        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
3267        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
3268        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
3269        buftubebot=bordtubebot+tubebufthick
3270        buftubetop=bordtubetop-tubebufthick
3271       endif
3272       ions=index(controlcard,"IONS").gt.0
3273       call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
3274       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
3275       write(iout,*) "R_CUT_ELE=",r_cut_ele
3276       call reada(controlcard,"R_CUT_MART",r_cut_mart,15.0d0)
3277       call reada(controlcard,"LAMBDA_MART",rlamb_mart,0.3d0)
3278       call reada(controlcard,"R_CUT_ANG",r_cut_ang,4.2d0)
3279       check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
3280       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
3281       call readi(controlcard,'SYM',symetr,1)
3282       write (iout,*) "DISTCHAINMAX",distchainmax
3283       write (iout,*) "delta",delta
3284       write (iout,*) "einicheck",einicheck
3285       write (iout,*) "rescale_mode",rescale_modeW
3286       call flush(iout)
3287       bxfile=index(controlcard,"BXFILE").gt.0
3288       cxfile=index(controlcard,"CXFILE").gt.0
3289       if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile) &
3290        bxfile=.true.
3291       histfile=index(controlcard,"HISTFILE").gt.0
3292       histout=index(controlcard,"HISTOUT").gt.0
3293       entfile=index(controlcard,"ENTFILE").gt.0
3294       zscfile=index(controlcard,"ZSCFILE").gt.0
3295       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
3296       write (iout,*) "with_dihed_constr ",with_dihed_constr
3297       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3298       return
3299       end subroutine read_general_data
3300 !------------------------------------------------------------------------------
3301       subroutine read_efree(*)
3302 !
3303 ! Read molecular data
3304 !
3305 !      implicit none
3306 !      include 'DIMENSIONS'
3307 !      include 'DIMENSIONS.ZSCOPT'
3308 !      include 'DIMENSIONS.COMPAR'
3309 !      include 'DIMENSIONS.FREE'
3310 !      include 'COMMON.IOUNITS'
3311 !      include 'COMMON.TIME1'
3312 !      include 'COMMON.SBRIDGE'
3313 !      include 'COMMON.CONTROL'
3314 !      include 'COMMON.CHAIN'
3315 !      include 'COMMON.HEADER'
3316 !      include 'COMMON.GEO'
3317 !      include 'COMMON.FREE'
3318       character(len=320) :: controlcard !,ucase
3319       integer :: iparm,ib,i,j,npars
3320 !el      integer ilen
3321 !el      external ilen
3322
3323       if (hamil_rep) then
3324         npars=1
3325       else
3326         npars=nParmSet
3327       endif
3328
3329 !      call alloc_wham_arrays
3330 !      allocate(nT_h(nParmSet))
3331 !      allocate(replica(nParmSet))
3332 !      allocate(umbrella(nParmSet))
3333 !      allocate(read_iset(nParmSet))
3334 !      allocate(nT_h(nParmSet))
3335
3336       do iparm=1,npars
3337
3338       call card_concat(controlcard,.true.)
3339       call readi(controlcard,'NT',nT_h(iparm),1)
3340       write (iout,*) "iparm",iparm," nt",nT_h(iparm)
3341       call flush(iout)
3342       if (nT_h(iparm).gt.MaxT_h) then
3343         write (iout,*)  "Error: parameter out of range: NT",nT_h(iparm),&
3344           MaxT_h
3345         return 1
3346       endif
3347       replica(iparm)=index(controlcard,"REPLICA").gt.0
3348       umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
3349       read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
3350       write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",&
3351         replica(iparm)," umbrella ",umbrella(iparm),&
3352         " read_iset",read_iset(iparm)
3353       call flush(iout)
3354       do ib=1,nT_h(iparm)
3355         call card_concat(controlcard,.true.)
3356         call readi(controlcard,'NR',nR(ib,iparm),1)
3357         if (umbrella(iparm)) then
3358           nRR(ib,iparm)=1
3359         else
3360           nRR(ib,iparm)=nR(ib,iparm)
3361         endif
3362         if (nR(ib,iparm).gt.MaxR) then
3363           write (iout,*)  "Error: parameter out of range: NR",&
3364             nR(ib,iparm),MaxR
3365         return 1
3366         endif
3367         call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
3368         beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
3369         call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),&
3370           0.0d0)
3371         do i=1,nR(ib,iparm)
3372           call card_concat(controlcard,.true.)
3373           call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,&
3374             100.0d0)
3375           call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,&
3376             0.0d0)
3377         enddo
3378       enddo
3379       do ib=1,nT_h(iparm)
3380         write (iout,*) "ib",ib," beta_h",&
3381           1.0d0/(0.001987*beta_h(ib,iparm))
3382         write (iout,*) "nR",nR(ib,iparm)
3383         write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
3384         do i=1,nR(ib,iparm)
3385           write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),&
3386             "q0",(q0(j,i,ib,iparm),j=1,nQ)
3387         enddo
3388         call flush(iout)
3389       enddo
3390
3391       enddo
3392
3393       if (hamil_rep) then
3394
3395        do iparm=2,nParmSet
3396           nT_h(iparm)=nT_h(1)
3397          do ib=1,nT_h(iparm)
3398            nR(ib,iparm)=nR(ib,1)
3399            if (umbrella(iparm)) then
3400              nRR(ib,iparm)=1
3401            else
3402              nRR(ib,iparm)=nR(ib,1)
3403            endif
3404            beta_h(ib,iparm)=beta_h(ib,1)
3405            do i=1,nR(ib,iparm)
3406              f(i,ib,iparm)=f(i,ib,1)
3407              do j=1,nQ
3408                KH(j,i,ib,iparm)=KH(j,i,ib,1)
3409                Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
3410              enddo
3411            enddo
3412            replica(iparm)=replica(1)
3413            umbrella(iparm)=umbrella(1)
3414            read_iset(iparm)=read_iset(1)
3415          enddo
3416        enddo
3417
3418       endif
3419
3420       return
3421       end subroutine read_efree
3422 !-----------------------------------------------------------------------------
3423       subroutine read_protein_data(*)
3424 !      implicit none
3425 !      include "DIMENSIONS"
3426 !      include "DIMENSIONS.ZSCOPT"
3427 !      include "DIMENSIONS.FREE"
3428 #ifdef MPI
3429       use MPI_data
3430       include "mpif.h"
3431       integer :: IERROR,ERRCODE!,STATUS(MPI_STATUS_SIZE)
3432 !      include "COMMON.MPI"
3433 #endif
3434 !      include "COMMON.CHAIN"
3435 !      include "COMMON.IOUNITS"
3436 !      include "COMMON.PROT"
3437 !      include "COMMON.PROTFILES"
3438 !      include "COMMON.NAMES"
3439 !      include "COMMON.FREE"
3440 !      include "COMMON.OBCINKA"
3441       character(len=64) :: nazwa
3442       character(len=16000) :: controlcard
3443       integer :: i,ii,ib,iR,iparm,nthr,npars !,ilen,iroof
3444 !el      external ilen,iroof
3445       if (hamil_rep) then
3446         npars=1
3447       else
3448         npars=nparmset
3449       endif
3450
3451       do iparm=1,npars
3452
3453 ! Read names of files with conformation data.
3454       if (replica(iparm)) then
3455         nthr = 1
3456       else
3457         nthr = nT_h(iparm)
3458       endif
3459       do ib=1,nthr
3460       do ii=1,nRR(ib,iparm)
3461       write (iout,*) "Parameter set",iparm," temperature",ib,&
3462        " window",ii
3463       call flush(iout)
3464       call card_concat(controlcard,.true.)
3465       write (iout,*) controlcard(:ilen(controlcard))
3466       call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
3467       call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
3468       call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
3469       call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
3470       call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),&
3471        maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
3472       call reada(controlcard,"TIME_START",&
3473         time_start_collect(ii,ib,iparm),0.0d0)
3474       call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),&
3475         1.0d10)
3476       write (iout,*) "rec_start",rec_start(ii,ib,iparm),&
3477        " rec_end",rec_end(ii,ib,iparm)
3478       write (iout,*) "time_start",time_start_collect(ii,ib,iparm),&
3479        " time_end",time_end_collect(ii,ib,iparm)
3480       call flush(iout)
3481       if (replica(iparm)) then
3482         call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
3483         write (iout,*) "Number of trajectories",totraj(ii,iparm)
3484         call flush(iout)
3485       endif
3486       if (nfile_bin(ii,ib,iparm).lt.2 &
3487           .and. nfile_asc(ii,ib,iparm).eq.0 &
3488           .and. nfile_cx(ii,ib,iparm).eq.0) then
3489         write (iout,*) "Error - no action specified!"
3490         return 1
3491       endif
3492       if (nfile_bin(ii,ib,iparm).gt.0) then
3493         call card_concat(controlcard,.false.)
3494         call split_string(controlcard,protfiles(1,1,ii,ib,iparm),&
3495          maxfile_prot,nfile_bin(ii,ib,iparm))
3496 #ifdef DEBUG
3497         write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
3498         write(iout,*) (protfiles(i,1,ii,ib,iparm),&
3499           i=1,nfile_bin(ii,ib,iparm))
3500 #endif
3501       endif
3502       if (nfile_asc(ii,ib,iparm).gt.0) then
3503         call card_concat(controlcard,.false.)
3504         call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
3505          maxfile_prot,nfile_asc(ii,ib,iparm))
3506 #ifdef DEBUG
3507         write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
3508         write(iout,*) (protfiles(i,2,ii,ib,iparm),&
3509           i=1,nfile_asc(ii,ib,iparm))
3510 #endif
3511       else if (nfile_cx(ii,ib,iparm).gt.0) then
3512         call card_concat(controlcard,.false.)
3513         call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
3514          maxfile_prot,nfile_cx(ii,ib,iparm))
3515 #ifdef DEBUG
3516         write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
3517         write(iout,*) (protfiles(i,2,ii,ib,iparm),&
3518          i=1,nfile_cx(ii,ib,iparm))
3519 #endif
3520       endif
3521       call flush(iout)
3522       enddo
3523       enddo
3524
3525       enddo
3526       return
3527       end subroutine read_protein_data
3528 !-------------------------------------------------------------------------------
3529       subroutine readsss(rekord,lancuch,wartosc,default)
3530 !      implicit none
3531       character*(*) :: rekord,lancuch,wartosc,default
3532       character(len=80) :: aux
3533       integer :: lenlan,lenrec,iread,ireade
3534 !el      external ilen
3535 !el      logical iblnk
3536 !el      external iblnk
3537       lenlan=ilen(lancuch)
3538       lenrec=ilen(rekord)
3539       iread=index(rekord,lancuch(:lenlan)//"=")
3540 !      print *,"rekord",rekord," lancuch",lancuch
3541 !      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
3542       if (iread.eq.0) then
3543         wartosc=default
3544         return
3545       endif
3546       iread=iread+lenlan+1
3547 !      print *,"iread",iread
3548 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3549       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
3550         iread=iread+1
3551 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3552       enddo
3553 !      print *,"iread",iread
3554       if (iread.gt.lenrec) then
3555          wartosc=default
3556         return
3557       endif
3558       ireade=iread+1
3559 !      print *,"ireade",ireade
3560       do while (ireade.lt.lenrec .and. &
3561          .not.iblnk(rekord(ireade:ireade)))
3562         ireade=ireade+1
3563       enddo
3564       wartosc=rekord(iread:ireade)
3565       return
3566       end subroutine readsss
3567 !----------------------------------------------------------------------------
3568       subroutine multreads(rekord,lancuch,tablica,dim,default)
3569 !      implicit none
3570       integer :: dim,i
3571       character*(*) rekord,lancuch,tablica(dim),default
3572       character(len=80) :: aux
3573       integer :: lenlan,lenrec,iread,ireade
3574 !el      external ilen
3575 !el      logical iblnk
3576 !el      external iblnk
3577       do i=1,dim
3578         tablica(i)=default
3579       enddo
3580       lenlan=ilen(lancuch)
3581       lenrec=ilen(rekord)
3582       iread=index(rekord,lancuch(:lenlan)//"=")
3583 !      print *,"rekord",rekord," lancuch",lancuch
3584 !      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
3585       if (iread.eq.0) return
3586       iread=iread+lenlan+1
3587       do i=1,dim
3588 !      print *,"iread",iread
3589 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3590       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
3591         iread=iread+1
3592 !      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
3593       enddo
3594 !      print *,"iread",iread
3595       if (iread.gt.lenrec) return
3596       ireade=iread+1
3597 !      print *,"ireade",ireade
3598       do while (ireade.lt.lenrec .and. &
3599          .not.iblnk(rekord(ireade:ireade)))
3600         ireade=ireade+1
3601       enddo
3602       tablica(i)=rekord(iread:ireade)
3603       iread=ireade+1
3604       enddo
3605       end subroutine multreads
3606 !----------------------------------------------------------------------------
3607       subroutine split_string(rekord,tablica,dim,nsub)
3608 !      implicit none
3609       integer :: dim,nsub,i,ii,ll,kk
3610       character*(*) tablica(dim)
3611       character*(*) rekord
3612 !el      integer ilen
3613 !el      external ilen
3614       do i=1,dim
3615         tablica(i)=" "
3616       enddo
3617       ii=1
3618       ll = ilen(rekord)
3619       nsub=0
3620       do i=1,dim
3621 ! Find the start of term name
3622         kk = 0
3623         do while (ii.le.ll .and. rekord(ii:ii).eq." ")
3624           ii = ii+1
3625         enddo
3626 ! Parse the name into TABLICA(i) until blank found
3627         do while (ii.le.ll .and. rekord(ii:ii).ne." ")
3628           kk = kk+1
3629           tablica(i)(kk:kk)=rekord(ii:ii)
3630           ii = ii+1
3631         enddo
3632         if (kk.gt.0) nsub=nsub+1
3633         if (ii.gt.ll) return
3634       enddo
3635       return
3636       end subroutine split_string
3637 !--------------------------------------------------------------------------------
3638 ! readrtns_compar.F
3639 !--------------------------------------------------------------------------------
3640       subroutine read_compar
3641 !
3642 ! Read molecular data
3643 !
3644       use conform_compar, only:alloc_compar_arrays
3645       use control_data, only:pdbref
3646       use geometry_data, only:deg2rad,rad2deg
3647 !      implicit none
3648 !      include 'DIMENSIONS'
3649 !      include 'DIMENSIONS.ZSCOPT'
3650 !      include 'DIMENSIONS.COMPAR'
3651 !      include 'DIMENSIONS.FREE'
3652 !      include 'COMMON.IOUNITS'
3653 !      include 'COMMON.TIME1'
3654 !      include 'COMMON.SBRIDGE'
3655 !      include 'COMMON.CONTROL'
3656 !      include 'COMMON.COMPAR'
3657 !      include 'COMMON.CHAIN'
3658 !      include 'COMMON.HEADER'
3659 !      include 'COMMON.GEO'
3660 !      include 'COMMON.FREE'
3661       character(len=320) :: controlcard !,ucase
3662       character(len=64) :: wfile
3663 !el      integer ilen
3664 !el      external ilen
3665       integer :: i,j,k
3666 !elwrite(iout,*)"jestesmy w read_compar"
3667       call card_concat(controlcard,.true.)
3668       pdbref=(index(controlcard,'PDBREF').gt.0)
3669       call reada(controlcard,'CUTOFF_UP',rmscut_base_up,4.0d0)
3670       call reada(controlcard,'CUTOFF_LOW',rmscut_base_low,3.0d0)
3671       call reada(controlcard,'RMSUP_LIM',rmsup_lim,4.0d0)
3672       call reada(controlcard,'RMSUPUP_LIM',rmsupup_lim,7.5d0)
3673       verbose = index(controlcard,"VERBOSE").gt.0
3674       lgrp=index(controlcard,"STATIN").gt.0
3675       lgrp_out=index(controlcard,"STATOUT").gt.0
3676       merge_helices=index(controlcard,"DONT_MERGE_HELICES").eq.0
3677       binary = index(controlcard,"BINARY").gt.0
3678       rmscut_base_up=rmscut_base_up/50
3679       rmscut_base_low=rmscut_base_low/50
3680       call reada(controlcard,"FRAC_SEC",frac_sec,0.66666666d0)
3681       call readi(controlcard,'NLEVEL',nlevel,1)
3682       if (nlevel.lt.0) then
3683         allocate(nfrag(2))
3684         call alloc_compar_arrays(maxfrag,1)
3685         goto 121
3686       else
3687         allocate(nfrag(nlevel))
3688       endif
3689 ! Read the data pertaining to elementary fragments (level 1)
3690       call readi(controlcard,'NFRAG',nfrag(1),0)
3691       write(iout,*)"nfrag(1)",nfrag(1)
3692       call alloc_compar_arrays(nfrag(1),nlevel)
3693       do j=1,nfrag(1)
3694         call card_concat(controlcard,.true.)
3695         write (iout,*) controlcard(:ilen(controlcard))
3696         call readi(controlcard,'NPIECE',npiece(j,1),0)
3697         call readi(controlcard,'N_SHIFT1',n_shift(1,j,1),0)
3698         call readi(controlcard,'N_SHIFT2',n_shift(2,j,1),0)
3699         call reada(controlcard,'ANGCUT',ang_cut(j),50.0d0)
3700         call reada(controlcard,'MAXANG',ang_cut1(j),360.0d0)
3701         call reada(controlcard,'FRAC_MIN',frac_min(j),0.666666d0)
3702         call reada(controlcard,'NC_FRAC',nc_fragm(j,1),0.5d0)
3703         call readi(controlcard,'NC_REQ',nc_req_setf(j,1),0)
3704         call readi(controlcard,'RMS',irms(j,1),0)
3705         call readi(controlcard,'LOCAL',iloc(j),1)
3706         call readi(controlcard,'ELCONT',ielecont(j,1),1)
3707         if (ielecont(j,1).eq.0) then
3708           call readi(controlcard,'SCCONT',isccont(j,1),1)
3709         endif
3710         ang_cut(j)=ang_cut(j)*deg2rad
3711         ang_cut1(j)=ang_cut1(j)*deg2rad
3712         do k=1,npiece(j,1)
3713           call card_concat(controlcard,.true.)
3714           call readi(controlcard,'IFRAG1',ifrag(1,k,j),0)
3715           call readi(controlcard,'IFRAG2',ifrag(2,k,j),0)
3716         enddo
3717         write(iout,*)"j",j," npiece",npiece(j,1)," ifrag",&
3718           (ifrag(1,k,j),ifrag(2,k,j),&
3719          k=1,npiece(j,1))," ang_cut",ang_cut(j)*rad2deg,&
3720           " ang_cut1",ang_cut1(j)*rad2deg
3721         write(iout,*)"n_shift",n_shift(1,j,1),n_shift(2,j,1)
3722         write(iout,*)"nc_frac",nc_fragm(j,1)," nc_req",nc_req_setf(j,1)
3723         write(iout,*)"irms",irms(j,1)," ielecont",ielecont(j,1),&
3724           " ilocal",iloc(j)," isccont",isccont(j,1)
3725       enddo
3726 ! Read data pertaning to higher levels
3727       do i=2,nlevel
3728         call card_concat(controlcard,.true.)
3729         call readi(controlcard,'NFRAG',NFRAG(i),0)
3730         write (iout,*) "i",i," nfrag",nfrag(i)
3731         do j=1,nfrag(i)
3732           call card_concat(controlcard,.true.)
3733           if (i.eq.2) then
3734             call readi(controlcard,'ELCONT',ielecont(j,i),0)
3735             if (ielecont(j,i).eq.0) then
3736               call readi(controlcard,'SCCONT',isccont(j,i),1)
3737             endif
3738             call readi(controlcard,'RMS',irms(j,i),0)
3739           else
3740             ielecont(j,i)=0
3741             isccont(j,i)=0
3742             irms(j,i)=1
3743           endif
3744           call readi(controlcard,'NPIECE',npiece(j,i),0)
3745           call readi(controlcard,'N_SHIFT1',n_shift(1,j,i),0)
3746           call readi(controlcard,'N_SHIFT2',n_shift(2,j,i),0)
3747           call multreadi(controlcard,'IPIECE',ipiece(1,j,i),&
3748             npiece(j,i),0)
3749           call reada(controlcard,'NC_FRAC',nc_fragm(j,i),0.5d0)
3750           call readi(controlcard,'NC_REQ',nc_req_setf(j,i),0)
3751           write(iout,*) "j",j," npiece",npiece(j,i)," n_shift",&
3752             n_shift(1,j,i),n_shift(2,j,i)," ielecont",ielecont(j,i),&
3753             " isccont",isccont(j,i)," irms",irms(j,i)
3754           write(iout,*) "ipiece",(ipiece(k,j,i),k=1,npiece(j,i))
3755           write(iout,*)"n_shift",n_shift(1,j,i),n_shift(2,j,i)
3756           write(iout,*)"nc_frac",nc_fragm(j,i),&
3757            " nc_req",nc_req_setf(j,i)
3758         enddo
3759       enddo
3760       if (binary) write (iout,*) "Classes written in binary format."
3761       return
3762   121 continue
3763       call reada(controlcard,'ANGCUT_HEL',angcut_hel,50.0d0)
3764       call reada(controlcard,'MAXANG_HEL',angcut1_hel,60.0d0)
3765       call reada(controlcard,'ANGCUT_BET',angcut_bet,90.0d0)
3766       call reada(controlcard,'MAXANG_BET',angcut1_bet,360.0d0)
3767       call reada(controlcard,'ANGCUT_STRAND',angcut_strand,90.0d0)
3768       call reada(controlcard,'MAXANG_STRAND',angcut1_strand,60.0d0)
3769       call reada(controlcard,'FRAC_MIN',frac_min_set,0.666666d0)
3770       call reada(controlcard,'NC_FRAC_HEL',ncfrac_hel,0.5d0)
3771       call readi(controlcard,'NC_REQ_HEL',ncreq_hel,0)
3772       call reada(controlcard,'NC_FRAC_BET',ncfrac_bet,0.5d0)
3773       call reada(controlcard,'NC_FRAC_PAIR',ncfrac_pair,0.3d0)
3774       call readi(controlcard,'NC_REQ_BET',ncreq_bet,0)
3775       call readi(controlcard,'NC_REQ_PAIR',ncreq_pair,0)
3776       call readi(controlcard,'NSHIFT_HEL',nshift_hel,3)
3777       call readi(controlcard,'NSHIFT_BET',nshift_bet,3)
3778       call readi(controlcard,'NSHIFT_STRAND',nshift_strand,3)
3779       call readi(controlcard,'NSHIFT_PAIR',nshift_pair,3)
3780       call readi(controlcard,'RMS_SINGLE',irms_single,0)
3781       call readi(controlcard,'CONT_SINGLE',icont_single,1)
3782       call readi(controlcard,'LOCAL_SINGLE',iloc_single,1)
3783       call readi(controlcard,'RMS_PAIR',irms_pair,0)
3784       call readi(controlcard,'CONT_PAIR',icont_pair,1)
3785       call readi(controlcard,'SPLIT_BET',isplit_bet,0)
3786       angcut_hel=angcut_hel*deg2rad
3787       angcut1_hel=angcut1_hel*deg2rad
3788       angcut_bet=angcut_bet*deg2rad
3789       angcut1_bet=angcut1_bet*deg2rad
3790       angcut_strand=angcut_strand*deg2rad
3791       angcut1_strand=angcut1_strand*deg2rad
3792       write (iout,*) "Automatic detection of structural elements"
3793       write (iout,*) 'NC_FRAC_HEL',ncfrac_hel,' NC_REQ_HEL',ncreq_hel,&
3794                      ' NC_FRAC_BET',ncfrac_bet,' NC_REQ_BET',ncreq_bet,&
3795                  ' RMS_SINGLE',irms_single,' CONT_SINGLE',icont_single,&
3796                  ' NC_FRAC_PAIR',ncfrac_pair,' NC_REQ_PAIR',ncreq_pair,&
3797         ' RMS_PAIR',irms_pair,' CONT_PAIR',icont_pair,&
3798         ' SPLIT_BET',isplit_bet
3799       write (iout,*) 'NSHIFT_HEL',nshift_hel,' NSHIFT_BET',nshift_bet,&
3800         ' NSHIFT_STRAND',nshift_strand,' NSHIFT_PAIR',nshift_pair
3801       write (iout,*) 'ANGCUT_HEL',angcut_hel*rad2deg,&
3802         ' MAXANG_HEL',angcut1_hel*rad2deg
3803       write (iout,*) 'ANGCUT_BET',angcut_bet*rad2deg,&
3804                      ' MAXANG_BET',angcut1_bet*rad2deg
3805       write (iout,*) 'ANGCUT_STRAND',angcut_strand*rad2deg,&
3806                      ' MAXANG_STRAND',angcut1_strand*rad2deg
3807       write (iout,*) 'FRAC_MIN',frac_min_set
3808       return
3809       end subroutine read_compar
3810 !--------------------------------------------------------------------------------
3811 ! read_ref_str.F
3812 !--------------------------------------------------------------------------------
3813       subroutine read_ref_structure(*)
3814 !
3815 ! Read the reference structure from the PDB file or from a PDB file or in the form of the dihedral
3816 ! angles.
3817 !
3818       use control_data, only:pdbref
3819       use geometry_data, only:nres,cref,c,dc,nsup,dc_norm,nend_sup,&
3820                               nstart_sup,nstart_seq,nperm,nres0
3821       use energy_data, only:nct,nnt,icont_ref,ncont_ref,itype,molnum
3822       use compare, only:seq_comp !,contact,elecont
3823       use geometry, only:chainbuild,dist
3824       use io_config, only:readpdb
3825 !
3826       use conform_compar, only:contact,elecont
3827 !      implicit none
3828 !      include 'DIMENSIONS'
3829 !      include 'DIMENSIONS.ZSCOPT'
3830 !      include 'DIMENSIONS.COMPAR'
3831 !      include 'COMMON.IOUNITS'
3832 !      include 'COMMON.GEO'
3833 !      include 'COMMON.VAR'
3834 !      include 'COMMON.INTERACT'
3835 !      include 'COMMON.LOCAL'
3836 !      include 'COMMON.NAMES'
3837 !      include 'COMMON.CHAIN'
3838 !      include 'COMMON.FFIELD'
3839 !      include 'COMMON.SBRIDGE'
3840 !      include 'COMMON.HEADER'
3841 !      include 'COMMON.CONTROL'
3842 !      include 'COMMON.CONTACTS1'
3843 !      include 'COMMON.PEPTCONT'
3844 !      include 'COMMON.TIME1'
3845 !      include 'COMMON.COMPAR'
3846       character(len=4) :: sequence(nres)
3847 !el      integer rescode
3848 !el      real(kind=8) :: x(maxvar)
3849       integer :: itype_pdb(nres,5)
3850 !el      logical seq_comp
3851       integer :: i,j,k,nres_pdb,iaux,mnum
3852       real(kind=8) :: ddsc !el,dist
3853       integer :: kkk !,ilen
3854 !el      external ilen
3855 !
3856       nres0=nres
3857       write (iout,*) "pdbref",pdbref
3858       if (pdbref) then
3859         read(inp,'(a)') pdbfile
3860         write (iout,'(2a,1h.)') 'PDB data will be read from file ',&
3861           pdbfile(:ilen(pdbfile))
3862         open(ipdbin,file=pdbfile,status='old',err=33)
3863         goto 34
3864   33    write (iout,'(a)') 'Error opening PDB file.'
3865         return 1
3866   34    continue
3867         do i=1,nres
3868           mnum=molnum(i)
3869           itype_pdb(i,mnum)=itype(i,mnum)
3870         enddo
3871
3872         call readpdb
3873
3874         do i=1,nres
3875           iaux=itype_pdb(i,mnum)
3876           itype_pdb(i,mnum)=itype(i,mnum)
3877           itype(i,mnum)=iaux
3878         enddo
3879         close (ipdbin)
3880         do kkk=1,nperm
3881         nres_pdb=nres
3882         nres=nres0
3883         nstart_seq=nnt
3884         if (nsup.le.(nct-nnt+1)) then
3885           do i=0,nct-nnt+1-nsup
3886             if (seq_comp(itype(nnt+i,molnum(nnt+i)),itype_pdb(nstart_sup,molnum(nstart_sup)),&
3887               nsup)) then
3888               do j=nnt+nsup-1,nnt,-1
3889                 do k=1,3
3890                   cref(k,nres+j+i,kkk)=cref(k,nres_pdb+j,kkk)
3891                 enddo
3892               enddo
3893               do j=nnt+nsup-1,nnt,-1
3894                 do k=1,3
3895                   cref(k,j+i,kkk)=cref(k,j,kkk)
3896                 enddo
3897                 write(iout,*) "J",j,"J+I",j+i
3898                 phi_ref(j+i)=phi_ref(j)
3899                 theta_ref(j+i)=theta_ref(j)
3900                 alph_ref(j+i)=alph_ref(j)
3901                 omeg_ref(j+i)=omeg_ref(j)
3902               enddo
3903 #ifdef DEBUG
3904               do j=nnt,nct
3905                 write (iout,'(i5,3f10.5,5x,3f10.5)') &
3906                   j,(cref(k,j,kkk),k=1,3),(cref(k,j+nres,kkk),k=1,3)
3907               enddo
3908 #endif
3909               nstart_seq=nnt+i
3910               nstart_sup=nnt+i
3911               goto 111
3912             endif
3913           enddo
3914           write (iout,'(a)') &
3915                   'Error - sequences to be superposed do not match.'
3916           return 1
3917         else
3918           do i=0,nsup-(nct-nnt+1)
3919             if (seq_comp(itype(nnt,molnum(nnt)),itype_pdb(nstart_sup+i,molnum(nstart_sup+i)),&
3920               nct-nnt+1)) &
3921             then
3922               nstart_sup=nstart_sup+i
3923               nsup=nct-nnt+1
3924               goto 111
3925             endif
3926           enddo
3927           write (iout,'(a)') &
3928                   'Error - sequences to be superposed do not match.'
3929         endif
3930         enddo
3931   111   continue
3932         write (iout,'(a,i5)') &
3933          'Experimental structure begins at residue',nstart_seq
3934       else
3935         call read_angles(inp,*38)
3936         goto 39
3937    38   write (iout,'(a)') 'Error reading reference structure.'
3938         return 1
3939    39   call chainbuild
3940         kkk=1
3941         nstart_sup=nnt
3942         nstart_seq=nnt
3943         nsup=nct-nnt+1
3944         do i=1,2*nres
3945           do j=1,3
3946             cref(j,i,kkk)=c(j,i)
3947           enddo
3948         enddo
3949       endif
3950       nend_sup=nstart_sup+nsup-1
3951       do i=1,2*nres
3952         do j=1,3
3953           c(j,i)=cref(j,i,kkk)
3954         enddo
3955       enddo
3956       do i=1,nres
3957         mnum=molnum(i)
3958         do j=1,3
3959           dc(j,nres+i)=cref(j,nres+i,kkk)-cref(j,i,kkk)
3960         enddo
3961         if (itype(i,mnum).ne.10) then
3962           ddsc = dist(i,nres+i)
3963           do j=1,3
3964             dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
3965           enddo
3966         else
3967           do j=1,3
3968             dc_norm(j,nres+i)=0.0d0
3969           enddo
3970         endif
3971 !        write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
3972 !         " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
3973 !         dc_norm(3,nres+i)**2
3974         do j=1,3
3975           dc(j,i)=c(j,i+1)-c(j,i)
3976         enddo
3977         ddsc = dist(i,i+1)
3978         do j=1,3
3979           dc_norm(j,i)=dc(j,i)/ddsc
3980         enddo
3981       enddo
3982 !      print *,"Calling contact"
3983       call contact(.true.,ncont_ref,icont_ref(1,1),&
3984         nstart_sup,nend_sup)
3985 !      print *,"Calling elecont"
3986       call elecont(.true.,ncont_pept_ref,&
3987          icont_pept_ref(1,1),&
3988          nstart_sup,nend_sup)
3989        write (iout,'(a,i3,a,i3,a,i3,a)') &
3990           'Number of residues to be superposed:',nsup,&
3991           ' (from residue',nstart_sup,' to residue',&
3992           nend_sup,').'
3993       return
3994       end subroutine read_ref_structure
3995 !--------------------------------------------------------------------------------
3996 ! geomout.F
3997 !--------------------------------------------------------------------------------
3998       subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev)
3999
4000       use geometry_data, only:nres,c,boxxsize,boxysize,boxzsize
4001       use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype,molnum
4002       use energy, only:boxshift
4003 !      implicit real*8 (a-h,o-z)
4004 !      include 'DIMENSIONS'
4005 !      include 'DIMENSIONS.ZSCOPT'
4006 !      include 'COMMON.CHAIN'
4007 !      include 'COMMON.INTERACT'
4008 !      include 'COMMON.NAMES'
4009 !      include 'COMMON.IOUNITS'
4010 !      include 'COMMON.HEADER'
4011 !      include 'COMMON.SBRIDGE'
4012       character(len=50) :: tytul
4013       character(len=1),dimension(24) :: chainid=reshape((/'A','B','C',&
4014                       'D','E','F','G','H','I','J','K','L','M','N','O',&
4015                       'P','Q','R','S','V','W','X','Y','Z'/),shape(chainid))
4016       integer,dimension(nres) :: ica !(maxres)
4017       real(kind=8) :: temp,efree,etot,entropy,rmsdev,xj,yj,zj
4018       integer :: ii,i,j,iti,ires,iatom,ichain,mnum
4019       write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')&
4020         ii,temp,rmsdev
4021       write (ipdb,'("REMARK DIMENSIONLESS FREE ENERGY",1pe15.5)') &
4022         efree
4023       write (ipdb,'("REMARK ENERGY",1pe15.5," ENTROPY",1pe15.5)') &
4024         etot,entropy
4025       iatom=0
4026       ichain=1
4027       ires=0
4028       do i=nnt,nct
4029         mnum=molnum(i)
4030         iti=itype(i,mnum)
4031         if (iti.eq.ntyp1_molec(mnum)) then
4032           if (itype(i-1,molnum(i-1)).eq.ntyp1_molec(mnum)) then
4033           ichain=ichain+1
4034           ires=0
4035           write (ipdb,'(a)') 'TER'
4036           endif
4037         else
4038         ires=ires+1
4039         iatom=iatom+1
4040         ica(i)=iatom
4041         if (mnum.ne.5) then
4042         write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
4043            ires,(c(j,i),j=1,3)
4044         else
4045         xj=boxshift(c(1,i)-c(1,2),boxxsize)
4046         yj=boxshift(c(2,i)-c(2,2),boxysize)
4047         zj=boxshift(c(3,i)-c(3,2),boxzsize)
4048         write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
4049            ires,c(1,2)+xj,c(2,2)+yj,c(3,2)+zj
4050         endif
4051         if ((iti.ne.10).and.(mnum.ne.5)) then
4052           iatom=iatom+1
4053           write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
4054             ires,(c(j,nres+i),j=1,3)
4055         endif
4056         endif
4057       enddo
4058       write (ipdb,'(a)') 'TER'
4059       do i=nnt,nct-1
4060         mnum=molnum(i)
4061         if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
4062         if (itype(i,mnum).eq.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
4063           write (ipdb,30) ica(i),ica(i+1)
4064         else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
4065           write (ipdb,30) ica(i),ica(i+1),ica(i)+1
4066         else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).eq.ntyp1_molec(mnum)) then
4067           write (ipdb,30) ica(i),ica(i)+1
4068         endif
4069       enddo
4070       if (itype(nct,molnum(nct)).ne.10) then
4071         write (ipdb,30) ica(nct),ica(nct)+1
4072       endif
4073       do i=1,nss
4074         write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
4075       enddo
4076       write (ipdb,'(a6)') 'ENDMDL'
4077   10  FORMAT ('ATOM',I7,'  CA  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
4078   20  FORMAT ('ATOM',I7,'  CB  ',A3,1X,A1,I4,4X,3F8.3,f15.3)
4079   30  FORMAT ('CONECT',8I5)
4080       return
4081       end subroutine pdboutW
4082 #endif
4083
4084       subroutine read_constr_homology
4085       use energy_data
4086       use control, only:init_int_table,homology_partition
4087       use MD_data, only:iset
4088       use geometry_data !only:nres,deg2rad,c,dc,nres_molec,crefjlee,cref
4089       use MPI_data, only:kolor
4090       use io_config, only:readpdb_template,readpdb
4091 !      implicit none
4092 !      include 'DIMENSIONS'
4093 !#ifdef MPI
4094 !      include 'mpif.h'
4095 !#endif
4096 !      include 'COMMON.SETUP'
4097 !      include 'COMMON.CONTROL'
4098 !      include 'COMMON.HOMOLOGY'
4099 !      include 'COMMON.CHAIN'
4100 !      include 'COMMON.IOUNITS'
4101 !      include 'COMMON.MD'
4102 !      include 'COMMON.QRESTR'
4103 !      include 'COMMON.GEO'
4104 !      include 'COMMON.INTERACT'
4105 !      include 'COMMON.NAMES'
4106 !      include 'COMMON.VAR'
4107 !
4108
4109 !     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
4110 !    &                 dist_cut
4111 !     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
4112 !    &    sigma_odl_temp(maxres,maxres,max_template)
4113       character*2 kic2
4114       character*24 model_ki_dist, model_ki_angle
4115       character*500 controlcard
4116       integer :: ki,i,ii,j,k,l
4117       integer, dimension (:), allocatable :: ii_in_use
4118       integer :: i_tmp,idomain_tmp,&
4119       irec,ik,iistart,nres_temp
4120 !      integer :: iset
4121 !      external :: ilen
4122       logical :: liiflag,lfirst
4123       integer :: i01,i10
4124 !
4125 !     FP - Nov. 2014 Temporary specifications for new vars
4126 !
4127       real(kind=8) :: rescore_tmp,x12,y12,z12,rescore2_tmp,&
4128                        rescore3_tmp, dist_cut
4129       real(kind=8), dimension (:,:),allocatable :: rescore
4130       real(kind=8), dimension (:,:),allocatable :: rescore2
4131       real(kind=8), dimension (:,:),allocatable :: rescore3
4132       real(kind=8) :: distal
4133       character*24 tpl_k_rescore
4134       character*256 pdbfile
4135
4136 ! -----------------------------------------------------------------
4137 ! Reading multiple PDB ref structures and calculation of retraints
4138 ! not using pre-computed ones stored in files model_ki_{dist,angle}
4139 ! FP (Nov., 2014)
4140 ! -----------------------------------------------------------------
4141 !
4142 !
4143 ! Alternative: reading from input
4144       call card_concat(controlcard,.true.)
4145       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
4146       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
4147       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
4148       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
4149       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
4150       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
4151       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
4152       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
4153       start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
4154 !      if(.not.read2sigma.and.start_from_model) then
4155 !          if(me1.eq.king .or. .not. out1file .and. fg_rank.eq.0)&
4156 !           write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
4157 !          start_from_model=.false.
4158 !          iranconf=(indpdb.le.0)
4159 !      endif
4160 !      if(start_from_model)&! .and. (me1.eq.king .or. .not. out1file))&
4161 !         write(iout,*) 'START_FROM_MODELS is ON'
4162 !      if(start_from_model .and. rest) then
4163 !        if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4164 !         write(iout,*) 'START_FROM_MODELS is OFF'
4165 !         write(iout,*) 'remove restart keyword from input'
4166 !        endif
4167 !      endif
4168 !      if (start_from_model) nmodel_start=constr_homology
4169       if(.not.allocated(waga_homology)) allocate (waga_homology(homol_nset))
4170       if (homol_nset.gt.1)then
4171          call card_concat(controlcard,.true.)
4172          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
4173 !         if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4174 !          write(iout,*) "iset homology_weight "
4175 !         do i=1,homol_nset
4176  !          write(iout,*) i,waga_homology(i)
4177  !         enddo
4178  !        endif
4179          iset=mod(kolor,homol_nset)+1
4180       else
4181        iset=1
4182        waga_homology(1)=1.0
4183       endif
4184
4185 !d      write (iout,*) "nnt",nnt," nct",nct
4186 !d      call flush(iout)
4187        if (.false.) then
4188        print *,"here klapaciuj"
4189 !      if (read_homol_frag) then
4190 !        call read_klapaucjusz
4191       else
4192
4193       lim_odl=0
4194       lim_dih=0
4195 !
4196 !      write(iout,*) 'nnt=',nnt,'nct=',nct
4197 !
4198 !      do i = nnt,nct
4199 !        do k=1,constr_homology
4200 !         idomain(k,i)=0
4201 !        enddo
4202 !      enddo
4203 !       idomain=0
4204
4205 !      ii=0
4206 !      do i = nnt,nct-2
4207 !        do j=i+2,nct
4208 !        ii=ii+1
4209 !        ii_in_use(ii)=0
4210 !        enddo
4211 !      enddo
4212       ii_in_use=0
4213       if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
4214       if(.not.allocated(chomo)) allocate(chomo(3,2*nres+2,constr_homology))
4215
4216       do k=1,constr_homology
4217
4218         read(inp,'(a)') pdbfile
4219         pdbfiles_chomo(k)=pdbfile
4220 !        if(me.eq.king .or. .not. out1file) &
4221          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',&
4222         pdbfile(:ilen(pdbfile))
4223         open(ipdbin,file=pdbfile,status='old',err=33)
4224         goto 34
4225   33    write (iout,'(a,5x,a)') 'Error opening PDB file',&
4226         pdbfile(:ilen(pdbfile))
4227         stop
4228   34    continue
4229 !        print *,'Begin reading pdb data'
4230 !
4231 ! Files containing res sim or local scores (former containing sigmas)
4232 !
4233
4234         write(kic2,'(bz,i2.2)') k
4235
4236         tpl_k_rescore="template"//kic2//".sco"
4237         write(iout,*) "tpl_k_rescore=",tpl_k_rescore
4238 !        unres_pdb=.false.
4239         nres_temp=nres
4240         write(iout,*) "read2sigma",read2sigma
4241
4242         if (read2sigma) then
4243           call readpdb_template(k)
4244         else
4245           call readpdb
4246         endif
4247         write(iout,*) "after readpdb"
4248         if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
4249         nres_chomo(k)=nres
4250         nres=nres_temp
4251         if(.not.allocated(rescore)) allocate(rescore(constr_homology,nres))
4252         if(.not.allocated(rescore2)) allocate(rescore2(constr_homology,nres))
4253         if(.not.allocated(rescore3)) allocate(rescore3(constr_homology,nres))
4254         if(.not.allocated(ii_in_use)) allocate(ii_in_use(nres*(nres-1)))
4255         if(.not.allocated(idomain)) allocate(idomain(constr_homology,nres))
4256         if(.not.allocated(l_homo)) allocate(l_homo(constr_homology,1000*nres))
4257         if(.not.allocated(ires_homo)) allocate(ires_homo(nres*200))
4258         if(.not.allocated(jres_homo)) allocate(jres_homo(nres*200))
4259         if(.not.allocated(odl)) allocate(odl(constr_homology,nres*200))
4260         if(.not.allocated(sigma_odl)) allocate(sigma_odl(constr_homology,nres*200))
4261         if(.not.allocated(dih)) allocate(dih(constr_homology,nres))
4262         if(.not.allocated(sigma_dih)) allocate(sigma_dih(constr_homology,nres))
4263         if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
4264         if(.not.allocated(sigma_theta)) allocate(sigma_theta(constr_homology,nres))
4265 !        if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
4266         if(.not.allocated(sigma_d)) allocate(sigma_d(constr_homology,nres))
4267         if(.not.allocated(xxtpl)) allocate(xxtpl(constr_homology,nres))
4268         if(.not.allocated(yytpl)) allocate(yytpl(constr_homology,nres))
4269         if(.not.allocated(zztpl)) allocate(zztpl(constr_homology,nres))
4270 !        if(.not.allocated(distance)) allocate(distance(constr_homology))
4271 !        if(.not.allocated(distancek)) allocate(distancek(constr_homology))
4272
4273
4274 !
4275 !     Distance restraints
4276 !
4277 !          ... --> odl(k,ii)
4278 ! Copy the coordinates from reference coordinates (?)
4279         do i=1,2*nres_chomo(k)
4280           do j=1,3
4281             c(j,i)=cref(j,i,1)
4282 !           write (iout,*) "c(",j,i,") =",c(j,i)
4283           enddo
4284         enddo
4285 !
4286 ! From read_dist_constr (commented out 25/11/2014 <-> res sim)
4287 !
4288 !         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
4289           open (ientin,file=tpl_k_rescore,status='old')
4290           if (nnt.gt.1) rescore(k,1)=0.0d0
4291           do irec=nnt,nct ! loop for reading res sim
4292             if (read2sigma) then
4293              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,&
4294                                      rescore3_tmp,idomain_tmp
4295              i_tmp=i_tmp+nnt-1
4296            write (*,*) "i_tmp", i_tmp,nnt
4297              idomain(k,i_tmp)=idomain_tmp
4298              rescore(k,i_tmp)=rescore_tmp
4299              rescore2(k,i_tmp)=rescore2_tmp
4300              rescore3(k,i_tmp)=rescore3_tmp
4301 !             if (.not. out1file .or. me.eq.king)&
4302 !             write(iout,'(a7,i5,3f10.5,i5)') "rescore",&
4303 !                           i_tmp,rescore2_tmp,rescore_tmp,&
4304 !                                     rescore3_tmp,idomain_tmp
4305             else
4306              idomain(k,irec)=1
4307              read (ientin,*,end=1401) rescore_tmp
4308
4309 !           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
4310              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
4311 !           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
4312             endif
4313           enddo
4314  1401   continue
4315         close (ientin)
4316         if (waga_dist.ne.0.0d0) then
4317           ii=0
4318           do i = nnt,nct-2
4319             do j=i+2,nct
4320
4321               x12=c(1,i)-c(1,j)
4322               y12=c(2,i)-c(2,j)
4323               z12=c(3,i)-c(3,j)
4324               distal=dsqrt(x12*x12+y12*y12+z12*z12)
4325 !              write (iout,*) k,i,j,distal,dist2_cut
4326
4327             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
4328                  .and. distal.le.dist2_cut ) then
4329
4330               ii=ii+1
4331               ii_in_use(ii)=1
4332               l_homo(k,ii)=.true.
4333
4334 !             write (iout,*) "k",k
4335 !             write (iout,*) "i",i," j",j," constr_homology",
4336 !    &                       constr_homology
4337               ires_homo(ii)=i
4338               jres_homo(ii)=j
4339               odl(k,ii)=distal
4340               if (read2sigma) then
4341                 sigma_odl(k,ii)=0
4342                 do ik=i,j
4343                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
4344                 enddo
4345                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
4346                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
4347               sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
4348               else
4349                 if (odl(k,ii).le.dist_cut) then
4350                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
4351                 else
4352 #ifdef OLDSIGMA
4353                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
4354                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
4355 #else
4356                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
4357                            dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
4358 #endif
4359                 endif
4360               endif
4361               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
4362             else
4363 !              ii=ii+1
4364 !              l_homo(k,ii)=.false.
4365             endif
4366             enddo
4367           enddo
4368         lim_odl=ii
4369         endif
4370 !        write (iout,*) "Distance restraints set"
4371 !        call flush(iout)
4372 !
4373 !     Theta, dihedral and SC retraints
4374 !
4375         if (waga_angle.gt.0.0d0) then
4376 !         open (ientin,file=tpl_k_sigma_dih,status='old')
4377 !         do irec=1,maxres-3 ! loop for reading sigma_dih
4378 !            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
4379 !            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
4380 !            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
4381 !    &                            sigma_dih(k,i+nnt-1)
4382 !         enddo
4383 !1402   continue
4384 !         close (ientin)
4385           do i = nnt+3,nct
4386             if (idomain(k,i).eq.0) then
4387                sigma_dih(k,i)=0.0
4388                cycle
4389             endif
4390             dih(k,i)=phiref(i) ! right?
4391 !           read (ientin,*) sigma_dih(k,i) ! original variant
4392 !             write (iout,*) "dih(",k,i,") =",dih(k,i)
4393 !             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
4394 !    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
4395 !    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
4396 !    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
4397
4398             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
4399                           rescore(k,i-2)+rescore(k,i-3))/4.0
4400 !            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
4401 !           write (iout,*) "Raw sigmas for dihedral angle restraints"
4402 !           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
4403 !           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
4404 !                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
4405 !   Instead of res sim other local measure of b/b str reliability possible
4406             if (sigma_dih(k,i).ne.0) &
4407             sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
4408 !           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
4409           enddo
4410           lim_dih=nct-nnt-2
4411         endif
4412 !        write (iout,*) "Dihedral angle restraints set"
4413 !        call flush(iout)
4414
4415         if (waga_theta.gt.0.0d0) then
4416 !         open (ientin,file=tpl_k_sigma_theta,status='old')
4417 !         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
4418 !            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
4419 !            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
4420 !    &                              sigma_theta(k,i+nnt-1)
4421 !         enddo
4422 !1403   continue
4423 !         close (ientin)
4424
4425           do i = nnt+2,nct ! right? without parallel.
4426 !         do i = i=1,nres ! alternative for bounds acc to readpdb?
4427 !         do i=ithet_start,ithet_end ! with FG parallel.
4428              if (idomain(k,i).eq.0) then
4429               sigma_theta(k,i)=0.0
4430               cycle
4431              endif
4432              thetatpl(k,i)=thetaref(i)
4433 !            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
4434 !            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
4435 !    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
4436 !    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
4437 !            read (ientin,*) sigma_theta(k,i) ! 1st variant
4438              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
4439                              rescore(k,i-2))/3.0
4440 !             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
4441              if (sigma_theta(k,i).ne.0) &
4442              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
4443
4444 !            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
4445 !                             rescore(k,i-2) !  right expression ?
4446 !            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
4447           enddo
4448         endif
4449 !        write (iout,*) "Angle restraints set"
4450 !        call flush(iout)
4451
4452         if (waga_d.gt.0.0d0) then
4453 !       open (ientin,file=tpl_k_sigma_d,status='old')
4454 !         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
4455 !            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
4456 !            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
4457 !    &                          sigma_d(k,i+nnt-1)
4458 !         enddo
4459 !1404   continue
4460
4461           do i = nnt,nct ! right? without parallel.
4462 !         do i=2,nres-1 ! alternative for bounds acc to readpdb?
4463 !         do i=loc_start,loc_end ! with FG parallel.
4464                if (itype(i,1).eq.10) cycle
4465                if (idomain(k,i).eq.0 ) then
4466                   sigma_d(k,i)=0.0
4467                   cycle
4468                endif
4469                xxtpl(k,i)=xxref(i)
4470                yytpl(k,i)=yyref(i)
4471                zztpl(k,i)=zzref(i)
4472 !              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
4473 !              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
4474 !              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
4475 !              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
4476                sigma_d(k,i)=rescore3(k,i) !  right expression ?
4477                if (sigma_d(k,i).ne.0) &
4478                sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
4479
4480 !              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
4481 !              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
4482 !              read (ientin,*) sigma_d(k,i) ! 1st variant
4483           enddo
4484         endif
4485       enddo
4486 !      write (iout,*) "SC restraints set"
4487 !      call flush(iout)
4488 !
4489 ! remove distance restraints not used in any model from the list
4490 ! shift data in all arrays
4491 !
4492 !      write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
4493       if (waga_dist.ne.0.0d0) then
4494         ii=0
4495         liiflag=.true.
4496         lfirst=.true.
4497         do i=nnt,nct-2
4498          do j=i+2,nct
4499           ii=ii+1
4500 !          if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
4501 !     &            .and. distal.le.dist2_cut ) then
4502 !          write (iout,*) "i",i," j",j," ii",ii
4503 !          call flush(iout)
4504           if (ii_in_use(ii).eq.0.and.liiflag.or. &
4505           ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
4506             liiflag=.false.
4507             i10=ii
4508             if (lfirst) then
4509               lfirst=.false.
4510               iistart=ii
4511             else
4512               if(i10.eq.lim_odl) i10=i10+1
4513               do ki=0,i10-i01-1
4514                ires_homo(iistart+ki)=ires_homo(ki+i01)
4515                jres_homo(iistart+ki)=jres_homo(ki+i01)
4516                ii_in_use(iistart+ki)=ii_in_use(ki+i01)
4517                do k=1,constr_homology
4518                 odl(k,iistart+ki)=odl(k,ki+i01)
4519                 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
4520                 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
4521                enddo
4522               enddo
4523               iistart=iistart+i10-i01
4524             endif
4525           endif
4526           if (ii_in_use(ii).ne.0.and..not.liiflag) then
4527              i01=ii
4528              liiflag=.true.
4529           endif
4530          enddo
4531         enddo
4532         lim_odl=iistart-1
4533       endif
4534 !      write (iout,*) "Removing distances completed"
4535 !      call flush(iout)
4536       endif ! .not. klapaucjusz
4537
4538       if (constr_homology.gt.0) call homology_partition
4539       write (iout,*) "After homology_partition"
4540       call flush(iout)
4541       if (constr_homology.gt.0) call init_int_table
4542       write (iout,*) "After init_int_table"
4543       call flush(iout)
4544 !      endif ! .not. klapaucjusz
4545 !      endif
4546 !      if (constr_homology.gt.0) call homology_partition
4547 !      write (iout,*) "After homology_partition"
4548 !      call flush(iout)
4549 !      if (constr_homology.gt.0) call init_int_table
4550 !      write (iout,*) "After init_int_table"
4551 !      call flush(iout)
4552 !      write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
4553 !      write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
4554 !
4555 ! Print restraints
4556 !
4557 #ifdef DEBUG
4558  !this debug needs correction
4559       if (.not.out_template_restr) return
4560 !d      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
4561       if(me1.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4562        write (iout,*) "Distance restraints from templates"
4563        do ii=1,lim_odl
4564        write(iout,'(3i7,100(2f8.2,1x,l1,4x))') &
4565         ii,ires_homo(ii),jres_homo(ii),&
4566         (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),&
4567         ki=1,constr_homology)
4568        enddo
4569        write (iout,*) "Dihedral angle restraints from templates"
4570        do i=nnt+3,nct
4571         write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
4572             (rad2deg*dih(ki,i),&
4573             rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
4574        enddo
4575        write (iout,*) "Virtual-bond angle restraints from templates"
4576        do i=nnt+2,nct
4577         write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
4578             (rad2deg*thetatpl(ki,i),&
4579             rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
4580        enddo
4581        write (iout,*) "SC restraints from templates"
4582        do i=nnt,nct
4583         write(iout,'(i7,100(4f8.2,4x))') i,&
4584         (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), &
4585          1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
4586        enddo
4587       endif
4588 #endif
4589       return
4590       end subroutine read_constr_homology
4591 !------------------------------------------------------------------------------
4592       end module io_wham
4593 !-----------------------------------------------------------------------------
4594 !-----------------------------------------------------------------------------
4595