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