source/wham/io_wham.F90

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