source/wham/io_wham.F90

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