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