source/wham/io_wham.F90

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