source/wham/io_wham.F90

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