source/wham/io_wham.F90

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