wham/io_wham.F90

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