source/unres/src-HCD-5D/readrtns_CSA.F

   1       subroutine readrtns
   2       implicit none
   3       include 'DIMENSIONS'
   4 #ifdef MPI
   5       include 'mpif.h'
   6 #endif
   7       include 'COMMON.SETUP'
   8       include 'COMMON.CONTROL'
   9       include 'COMMON.SBRIDGE'
  10       include 'COMMON.IOUNITS'
  11       include 'COMMON.SPLITELE'
  12       integer i,j
  13       logical file_exist
  14 C Read job setup parameters
  15       call read_control
  16 C Read force-field parameters except weights
  17       call parmread
  18 C Read control parameters for energy minimzation if required
  19       if (minim) call read_minim
  20 C Read MCM control parameters if required
  21       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
  22 C Read MD control parameters if reqjuired
  23       if (modecalc.eq.12) call read_MDpar
  24 C Read MREMD control parameters if required
  25       if (modecalc.eq.14) then
  26          call read_MDpar
  27          call read_REMDpar
  28       endif
  29 C Read MUCA control parameters if required
  30       if (lmuca) call read_muca
  31 C Read CSA control parameters if required (from fort.40 if exists
  32 C otherwise from general input file)
  33       if (modecalc.eq.8) then
  34        inquire (file="fort.40",exist=file_exist)
  35        if (.not.file_exist) call csaread
  36       endif
  37 cfmc      if (modecalc.eq.10) call mcmfread
  38 c Read energy-term weights and disulfide parameters
  39       call weightread
  40 C Read molecule information, molecule geometry, energy-term weights, and
  41 C restraints if requested
  42       call molread
  43 C Print restraint information
  44 #ifdef MPI
  45       if (.not. out1file .or. me.eq.king) then
  46 #endif
  47       if (nhpb.gt.nss)  then
  48         write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
  49      &  "The following",nhpb-nss,
  50      &  " distance restraints have been imposed:",
  51      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
  52      &  "  score"," type"
  53         do i=nss+1,nhpb
  54           write (iout,'(4i5,2f8.2,3f10.5,2i5)')i-nss,ihpb(i),jhpb(i),
  55      &  ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
  56      &  irestr_type(i)
  57         enddo
  58       endif
  59       if (npeak.gt.0) then
  60         write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
  61      &  "The following",npeak,
  62      &  " NMR peak restraints have been imposed:",
  63      &  "   Nr"," res1"," res2"," beta","   d1","   d2","    k","    V",
  64      &  "  score"," type"," ipeak"
  65         do i=1,npeak
  66           do j=ipeak(1,i),ipeak(2,i)
  67             write (iout,'(5i5,2f8.2,2f10.5,i5)')i,j,ihpb_peak(j),
  68      &       jhpb_peak(j),ibecarb_peak(j),dhpb_peak(j),dhpb1_peak(j),
  69      &       forcon_peak(j),fordepth_peak(i),irestr_type_peak(j)
  70           enddo
  71         enddo
  72         write (iout,*) "The ipeak array"
  73         do i=1,npeak
  74           write (iout,'(3i5)' ) i,ipeak(1,i),ipeak(2,i)
  75         enddo
  76       endif
  77 #ifdef MPI
  78       endif
  79 #endif
  80 c      print *,"Processor",myrank," leaves READRTNS"
  81       return
  82       end
  83 C-------------------------------------------------------------------------------
  84       subroutine read_control
  85 C
  86 C Read contorl data
  87 C
  88       implicit none
  89       include 'DIMENSIONS'
  90 #ifdef MP
  91       include 'mpif.h'
  92       logical OKRandom, prng_restart
  93       double precision  r1
  94 #endif
  95       include 'COMMON.IOUNITS'
  96       include 'COMMON.TIME1'
  97       include 'COMMON.THREAD'
  98       include 'COMMON.SBRIDGE'
  99       include 'COMMON.CONTROL'
 100       include 'COMMON.SAXS'
 101       include 'COMMON.MCM'
 102       include 'COMMON.MAP'
 103       include 'COMMON.HEADER'
 104       include 'COMMON.CSA'
 105       include 'COMMON.CHAIN'
 106       include 'COMMON.MUCA'
 107       include 'COMMON.MD'
 108       include 'COMMON.FFIELD'
 109       include 'COMMON.INTERACT'
 110       include 'COMMON.SETUP'
 111       include 'COMMON.SPLITELE'
 112       include 'COMMON.SHIELD'
 113       include 'COMMON.GEO'
 114       integer i
 115       integer KDIAG,ICORFL,IXDR
 116       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
 117       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
 118       character*80 ucase
 119       character*320 controlcard
 120       double precision seed
 121
 122       nglob_csa=0
 123       eglob_csa=1d99
 124       nmin_csa=0
 125       read (INP,'(a)') titel
 126       call card_concat(controlcard)
 127 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
 128 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
 129       call reada(controlcard,'SEED',seed,0.0D0)
 130       call random_init(seed)
 131 C Set up the time limit (caution! The time must be input in minutes!)
 132       read_cart=index(controlcard,'READ_CART').gt.0
 133       out_cart=index(controlcard,'OUT_CART').gt.0
 134       out_int=index(controlcard,'OUT_INT').gt.0
 135       gmatout=index(controlcard,'GMATOUT').gt.0
 136       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
 137 C this variable with_theta_constr is the variable which allow to read and execute the
 138 C constrains on theta angles WITH_THETA_CONSTR is the keyword
 139       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
 140       write (iout,*) "with_theta_constr ",with_theta_constr
 141       call readi(controlcard,'NSAXS',nsaxs,0)
 142       call readi(controlcard,'SAXS_MODE',saxs_mode,0)
 143       call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
 144       call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
 145       write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
 146      &   SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
 147       call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
 148       read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
 149       out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
 150       out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
 151       call readi(controlcard,'SYM',symetr,1)
 152       call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
 153       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
 154       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
 155 c      call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
 156 c      call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
 157 c      call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
 158 c      call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 159 c      call reada(controlcard,'DRMS',drms,0.1D0)
 160 c      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 161 c       write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 162 c       write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 163 c       write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 164 c       write (iout,'(a,f10.1)')'DRMS    = ',drms
 165 cc       write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 166 c       write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 167 c      endif
 168       call readi(controlcard,'NZ_START',nz_start,0)
 169       call readi(controlcard,'NZ_END',nz_end,0)
 170 c      call readi(controlcard,'IZ_SC',iz_sc,0)
 171       timlim=60.0D0*timlim
 172       safety = 60.0d0*safety
 173       modecalc=0
 174       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 175       minim=(index(controlcard,'MINIMIZE').gt.0)
 176       dccart=(index(controlcard,'CART').gt.0)
 177       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 178       overlapsc=.not.overlapsc
 179       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
 180       searchsc=.not.searchsc
 181       sideadd=(index(controlcard,'SIDEADD').gt.0)
 182       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 183       mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
 184       outpdb=(index(controlcard,'PDBOUT').gt.0)
 185       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 186       pdbref=(index(controlcard,'PDBREF').gt.0)
 187       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 188       indpdb=index(controlcard,'PDBSTART')
 189       extconf=(index(controlcard,'EXTCONF').gt.0)
 190       AFMlog=(index(controlcard,'AFM'))
 191       selfguide=(index(controlcard,'SELFGUIDE'))
 192 c      print *,'AFMlog',AFMlog,selfguide,"KUPA"
 193       call readi(controlcard,'TUBEMOD',tubelog,0)
 194 c      write (iout,*) TUBElog,"TUBEMODE"
 195       call readi(controlcard,'IPRINT',iprint,0)
 196 C SHIELD keyword sets if the shielding effect of side-chains is used
 197 C 0 denots no shielding is used all peptide are equally despite the
 198 C solvent accesible area
 199 C 1 the newly introduced function
 200 C 2 reseved for further possible developement
 201       call readi(controlcard,'SHIELD',shield_mode,0)
 202 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 203         write(iout,*) "shield_mode",shield_mode
 204 C      endif
 205       call readi(controlcard,'TORMODE',tor_mode,0)
 206 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 207         write(iout,*) "torsional and valence angle mode",tor_mode
 208       call readi(controlcard,'MAXGEN',maxgen,10000)
 209       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 210       call readi(controlcard,"KDIAG",kdiag,0)
 211       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
 212       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 213      & write (iout,*) "RESCALE_MODE",rescale_mode
 214       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 215       if (index(controlcard,'REGULAR').gt.0.0D0) then
 216         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 217         modecalc=1
 218         refstr=.true.
 219       endif
 220       if (refstr) then
 221       if (index(controlcard,"CASC").gt.0) then
 222         iz_sc=1
 223 c      else if (index(controlcard,"CAONLY").gt.0) then
 224 c        iz_sc=0
 225       else if (index(controlcard,"SCONLY").gt.0) then
 226         iz_sc=2
 227       else
 228         iz_sc=0
 229 c        write (iout,*) "ERROR! Must specify CASC CAONLY or SCONLY when",
 230 c     &    " specifying REFSTR, PDBREF or REGULAR."
 231 c        stop
 232       endif
 233       endif
 234       if (index(controlcard,'CHECKGRAD').gt.0) then
 235         modecalc=5
 236         if (index(controlcard,' CART').gt.0) then
 237           icheckgrad=1
 238         elseif (index(controlcard,'CARINT').gt.0) then
 239           icheckgrad=2
 240         else
 241           icheckgrad=3
 242         endif
 243         call reada(controlcard,'DELTA',aincr,1.0d-7)
 244 c        write (iout,*) "icheckgrad",icheckgrad
 245       elseif (index(controlcard,'THREAD').gt.0) then
 246         modecalc=2
 247         call readi(controlcard,'THREAD',nthread,0)
 248         if (nthread.gt.0) then
 249           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 250         else
 251           if (fg_rank.eq.0)
 252      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 253           stop 'Error termination in Read_Control.'
 254         endif
 255       else if (index(controlcard,'MCMA').gt.0) then
 256         modecalc=3
 257       else if (index(controlcard,'MCEE').gt.0) then
 258         modecalc=6
 259       else if (index(controlcard,'MULTCONF').gt.0) then
 260         modecalc=4
 261       else if (index(controlcard,'MAP').gt.0) then
 262         modecalc=7
 263         call readi(controlcard,'MAP',nmap,0)
 264       else if (index(controlcard,'CSA').gt.0) then
 265         modecalc=8
 266 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 267 crc
 268 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 269 crc
 270 crc        modecalc=9
 271 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 272 cfmc        modecalc=10
 273       else if (index(controlcard,'SOFTREG').gt.0) then
 274         modecalc=11
 275       else if (index(controlcard,'CHECK_BOND').gt.0) then
 276         modecalc=-1
 277       else if (index(controlcard,'TEST').gt.0) then
 278         modecalc=-2
 279       else if (index(controlcard,'MD').gt.0) then
 280         modecalc=12
 281       else if (index(controlcard,'RE ').gt.0) then
 282         modecalc=14
 283       endif
 284
 285       lmuca=index(controlcard,'MUCA').gt.0
 286       call readi(controlcard,'MUCADYN',mucadyn,0)
 287       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 288       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 289      & then
 290        write (iout,*) 'MUCADYN=',mucadyn
 291        write (iout,*) 'MUCASMOOTH=',muca_smooth
 292       endif
 293
 294       iscode=index(controlcard,'ONE_LETTER')
 295       indphi=index(controlcard,'PHI')
 296       indback=index(controlcard,'BACK')
 297       iranconf=index(controlcard,'RAND_CONF')
 298       i2ndstr=index(controlcard,'USE_SEC_PRED')
 299       gradout=index(controlcard,'GRADOUT').gt.0
 300       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 301 C  DISTCHAINMAX become obsolete for periodic boundry condition
 302       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
 303 C Reading the dimensions of box in x,y,z coordinates
 304       call reada(controlcard,'BOXX',boxxsize,100.0d0)
 305       call reada(controlcard,'BOXY',boxysize,100.0d0)
 306       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
 307 c Cutoff range for interactions
 308       call reada(controlcard,"R_CUT",r_cut,15.0d0)
 309       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 310       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
 311       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
 312       if (lipthick.gt.0.0d0) then
 313        bordliptop=(boxzsize+lipthick)/2.0
 314        bordlipbot=bordliptop-lipthick
 315 C      endif
 316       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
 317      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
 318       buflipbot=bordlipbot+lipbufthick
 319       bufliptop=bordliptop-lipbufthick
 320       if ((lipbufthick*2.0d0).gt.lipthick)
 321      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
 322       endif
 323       write(iout,*) "bordliptop=",bordliptop
 324       write(iout,*) "bordlipbot=",bordlipbot
 325       write(iout,*) "bufliptop=",bufliptop
 326       write(iout,*) "buflipbot=",buflipbot
 327       write (iout,*) "SHIELD MODE",shield_mode
 328       if (TUBElog.gt.0) then
 329        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
 330        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
 331        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
 332        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
 333        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
 334        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
 335        buftubebot=bordtubebot+tubebufthick
 336        buftubetop=bordtubetop-tubebufthick
 337       endif
 338       if (me.eq.king .or. .not.out1file )
 339      &  write (iout,*) "DISTCHAINMAX",distchainmax
 340
 341       if(me.eq.king.or..not.out1file)
 342      & write (iout,'(2a)') diagmeth(kdiag),
 343      &  ' routine used to diagonalize matrices.'
 344       return
 345       end
 346 c--------------------------------------------------------------------------
 347       subroutine read_REMDpar
 348 C
 349 C Read REMD settings
 350 C
 351       implicit none
 352       include 'DIMENSIONS'
 353       include 'COMMON.IOUNITS'
 354       include 'COMMON.TIME1'
 355       include 'COMMON.MD'
 356 #ifndef LANG0
 357       include 'COMMON.LANGEVIN'
 358 #else
 359 #ifdef FIVEDIAG
 360       include 'COMMON.LANGEVIN.lang0.5diag'
 361 #else
 362       include 'COMMON.LANGEVIN.lang0'
 363 #endif
 364 #endif
 365       include 'COMMON.INTERACT'
 366       include 'COMMON.NAMES'
 367       include 'COMMON.GEO'
 368       include 'COMMON.REMD'
 369       include 'COMMON.CONTROL'
 370       include 'COMMON.SETUP'
 371       character*80 ucase
 372       character*320 controlcard
 373       character*3200 controlcard1
 374       integer iremd_m_total,i
 375
 376       if(me.eq.king.or..not.out1file)
 377      & write (iout,*) "REMD setup"
 378
 379       call card_concat(controlcard)
 380       call readi(controlcard,"NREP",nrep,3)
 381       call readi(controlcard,"NSTEX",nstex,1000)
 382       call reada(controlcard,"RETMIN",retmin,10.0d0)
 383       call reada(controlcard,"RETMAX",retmax,1000.0d0)
 384       mremdsync=(index(controlcard,'SYNC').gt.0)
 385       call readi(controlcard,"NSYN",i_sync_step,100)
 386       restart1file=(index(controlcard,'REST1FILE').gt.0)
 387       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
 388       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
 389       if(max_cache_traj_use.gt.max_cache_traj)
 390      &           max_cache_traj_use=max_cache_traj
 391       if(me.eq.king.or..not.out1file) then
 392 cd       if (traj1file) then
 393 crc caching is in testing - NTWX is not ignored
 394 cd        write (iout,*) "NTWX value is ignored"
 395 cd        write (iout,*) "  trajectory is stored to one file by master"
 396 cd        write (iout,*) "  before exchange at NSTEX intervals"
 397 cd       endif
 398        write (iout,*) "NREP= ",nrep
 399        write (iout,*) "NSTEX= ",nstex
 400        write (iout,*) "SYNC= ",mremdsync
 401        write (iout,*) "NSYN= ",i_sync_step
 402        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
 403       endif
 404       remd_tlist=.false.
 405       if (index(controlcard,'TLIST').gt.0) then
 406          remd_tlist=.true.
 407          call card_concat(controlcard1)
 408          read(controlcard1,*) (remd_t(i),i=1,nrep)
 409          if(me.eq.king.or..not.out1file)
 410      &    write (iout,*)'tlist',(remd_t(i),i=1,nrep)
 411       endif
 412       remd_mlist=.false.
 413       if (index(controlcard,'MLIST').gt.0) then
 414          remd_mlist=.true.
 415          call card_concat(controlcard1)
 416          read(controlcard1,*) (remd_m(i),i=1,nrep)
 417          if(me.eq.king.or..not.out1file) then
 418           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
 419           iremd_m_total=0
 420           do i=1,nrep
 421            iremd_m_total=iremd_m_total+remd_m(i)
 422           enddo
 423            write (iout,*) 'Total number of replicas ',iremd_m_total
 424           endif
 425          endif
 426         if (adaptive) then
 427           write (iout,*)
 428      &      "Adaptive (PMF-biased) umbrella sampling will be run"
 429           call PMFread
 430         endif
 431       if(me.eq.king.or..not.out1file)
 432      &   write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
 433       return
 434       end
 435 c--------------------------------------------------------------------------
 436       subroutine read_MDpar
 437 C
 438 C Read MD settings
 439 C
 440       implicit none
 441       include 'DIMENSIONS'
 442       include 'COMMON.IOUNITS'
 443       include 'COMMON.TIME1'
 444       include 'COMMON.MD'
 445       include 'COMMON.QRESTR'
 446 #ifndef LANG0
 447       include 'COMMON.LANGEVIN'
 448 #else
 449 #ifdef FIVEDIAG
 450       include 'COMMON.LANGEVIN.lang0.5diag'
 451 #else
 452       include 'COMMON.LANGEVIN.lang0'
 453 #endif
 454 #endif
 455       include 'COMMON.INTERACT'
 456       include 'COMMON.NAMES'
 457       include 'COMMON.GEO'
 458       include 'COMMON.SETUP'
 459       include 'COMMON.CONTROL'
 460       include 'COMMON.SPLITELE'
 461       include 'COMMON.FFIELD'
 462       character*80 ucase
 463       character*320 controlcard
 464       integer i
 465       double precision eta
 466
 467       call card_concat(controlcard)
 468       call readi(controlcard,"NSTEP",n_timestep,1000000)
 469       call readi(controlcard,"NTWE",ntwe,100)
 470       call readi(controlcard,"NTWX",ntwx,1000)
 471       call reada(controlcard,"DT",d_time,1.0d-1)
 472       call reada(controlcard,"DVMAX",dvmax,2.0d1)
 473       call reada(controlcard,"DAMAX",damax,1.0d1)
 474       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
 475       call readi(controlcard,"LANG",lang,0)
 476       RESPA = index(controlcard,"RESPA") .gt. 0
 477       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
 478       ntime_split0=ntime_split
 479       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
 480       ntime_split0=ntime_split
 481 c      call reada(controlcard,"R_CUT",r_cut,2.0d0)
 482 c      call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 483       rest = index(controlcard,"REST").gt.0
 484       tbf = index(controlcard,"TBF").gt.0
 485       usampl = index(controlcard,"USAMPL").gt.0
 486       scale_umb =  index(controlcard,"SCALE_UMB").gt.0
 487       adaptive = index(controlcard,"ADAPTIVE").gt.0
 488       mdpdb = index(controlcard,"MDPDB").gt.0
 489       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 490       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
 491       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
 492       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
 493       if (count_reset_moment.eq.0) count_reset_moment=1000000000
 494       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
 495       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
 496       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
 497       if (count_reset_vel.eq.0) count_reset_vel=1000000000
 498       large = index(controlcard,"LARGE").gt.0
 499       print_compon = index(controlcard,"PRINT_COMPON").gt.0
 500       rattle = index(controlcard,"RATTLE").gt.0
 501       preminim = index(controlcard,"PREMINIM").gt.0
 502       if (iranconf.gt.0 .or. indpdb.gt.0 .or. start_from_model) then
 503         preminim=.true.
 504         dccart=.true.
 505       endif
 506       write (iout,*) "PREMINIM ",preminim
 507       if (preminim) then
 508         if (index(controlcard,'CART').gt.0) dccart=.true.
 509         write (iout,*) "dccart ",dccart
 510         write (iout,*) "read_minim ",index(controlcard,'READ_MINIM_PAR')
 511         if (index(controlcard,'READ_MINIM_PAR').gt.0) call read_minim
 512       endif
 513 c  if performing umbrella sampling, fragments constrained are read from the fragment file
 514       nset=0
 515       if(usampl) then
 516         write (iout,*) "Umbrella sampling will be run"
 517         if (scale_umb.and.adaptive) then
 518           write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
 519           write (iout,*) "Select one of those and re-run the job."
 520           stop
 521         endif
 522         if (scale_umb) write (iout,*)
 523      &"Umbrella-restraint force constants will be scaled by temperature"
 524         call read_fragments
 525       endif
 526
 527       if(me.eq.king.or..not.out1file) then
 528        write (iout,*)
 529        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
 530        write (iout,*)
 531        write (iout,'(a)') "The units are:"
 532        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
 533        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
 534      &  " acceleration: angstrom/(48.9 fs)**2"
 535        write (iout,'(a)') "energy: kcal/mol, temperature: K"
 536        write (iout,*)
 537        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
 538        write (iout,'(a60,f10.5,a)')
 539      &  "Initial time step of numerical integration:",d_time,
 540      &  " natural units"
 541        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
 542        if (RESPA) then
 543         write (iout,'(2a,i4,a)')
 544      &    "A-MTS algorithm used; initial time step for fast-varying",
 545      &    " short-range forces split into",ntime_split," steps."
 546         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
 547      &   r_cut," lambda",rlamb
 548        endif
 549        write (iout,'(2a,f10.5)')
 550      &  "Maximum acceleration threshold to reduce the time step",
 551      &  "/increase split number:",damax
 552        write (iout,'(2a,f10.5)')
 553      &  "Maximum predicted energy drift to reduce the timestep",
 554      &  "/increase split number:",edriftmax
 555        write (iout,'(a60,f10.5)')
 556      & "Maximum velocity threshold to reduce velocities:",dvmax
 557        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
 558        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
 559        if (rattle) write (iout,'(a60)')
 560      &  "Rattle algorithm used to constrain the virtual bonds"
 561       endif
 562       reset_fricmat=1000
 563       if (lang.gt.0) then
 564         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
 565         call reada(controlcard,"RWAT",rwat,1.4d0)
 566         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
 567         surfarea=index(controlcard,"SURFAREA").gt.0
 568         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
 569         if(me.eq.king.or..not.out1file)then
 570          write (iout,'(/a,$)') "Langevin dynamics calculation"
 571          if (lang.eq.1) then
 572           write (iout,'(a/)')
 573      &      " with direct integration of Langevin equations"
 574          else if (lang.eq.2) then
 575           write (iout,'(a/)') " with TINKER stochasic MD integrator"
 576          else if (lang.eq.3) then
 577           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
 578          else if (lang.eq.4) then
 579           write (iout,'(a/)') " in overdamped mode"
 580          else
 581           write (iout,'(//a,i5)')
 582      &      "=========== ERROR: Unknown Langevin dynamics mode:",lang
 583           stop
 584          endif
 585          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 586          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
 587          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
 588          write (iout,'(a60,f10.5)')
 589      &   "Scaling factor of the friction forces:",scal_fric
 590          if (surfarea) write (iout,'(2a,i10,a)')
 591      &     "Friction coefficients will be scaled by solvent-accessible",
 592      &     " surface area every",reset_fricmat," steps."
 593         endif
 594 c Calculate friction coefficients and bounds of stochastic forces
 595         eta=6*pi*cPoise*etawat
 596         if(me.eq.king.or..not.out1file)
 597      &   write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
 598      &   ,eta
 599         gamp=scal_fric*(pstok+rwat)*eta
 600         stdfp=dsqrt(2*Rb*t_bath/d_time)
 601         do i=1,ntyp
 602           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
 603           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 604         enddo
 605         if(me.eq.king.or..not.out1file)then
 606          write (iout,'(/2a/)')
 607      &   "Radii of site types and friction coefficients and std's of",
 608      &   " stochastic forces of fully exposed sites"
 609          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
 610          do i=1,ntyp
 611           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
 612      &     gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
 613          enddo
 614         endif
 615       else if (tbf) then
 616         if(me.eq.king.or..not.out1file)then
 617          write (iout,'(a)') "Berendsen bath calculation"
 618          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 619          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
 620          if (reset_moment)
 621      &   write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
 622      &   count_reset_moment," steps"
 623          if (reset_vel)
 624      &    write (iout,'(a,i10,a)')
 625      &    "Velocities will be reset at random every",count_reset_vel,
 626      &   " steps"
 627         endif
 628       else
 629         if(me.eq.king.or..not.out1file)
 630      &   write (iout,'(a31)') "Microcanonical mode calculation"
 631       endif
 632       if(me.eq.king.or..not.out1file)then
 633        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
 634        if (usampl) then
 635           write(iout,*) "MD running with constraints."
 636           write(iout,*) "Equilibration time ", eq_time, " mtus."
 637           write(iout,*) "Constraining ", nfrag," fragments."
 638           write(iout,*) "Length of each fragment, weight and q0:"
 639           do iset=1,nset
 640            write (iout,*) "Set of restraints #",iset
 641            do i=1,nfrag
 642               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
 643      &           ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
 644            enddo
 645            write(iout,*) "constraints between ", npair, "fragments."
 646            write(iout,*) "constraint pairs, weights and q0:"
 647            do i=1,npair
 648             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
 649      &             ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
 650            enddo
 651            write(iout,*) "angle constraints within ", nfrag_back,
 652      &      "backbone fragments."
 653            if (loc_qlike) then
 654            write(iout,*) "fragment, weights, q0:"
 655            do i=1,nfrag_back
 656             write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
 657      &         ifrag_back(2,i,iset),
 658      &         wfrag_back(1,i,iset),qin_back(1,i,iset),
 659      &         wfrag_back(2,i,iset),qin_back(2,i,iset),
 660      &         wfrag_back(3,i,iset),qin_back(3,i,iset)
 661            enddo
 662            else
 663            write(iout,*) "fragment, weights:"
 664            do i=1,nfrag_back
 665             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
 666      &         ifrag_back(2,i,iset),wfrag_back(1,i,iset),
 667      &         wfrag_back(2,i,iset),wfrag_back(3,i,iset)
 668            enddo
 669            endif
 670           enddo
 671         iset=mod(kolor,nset)+1
 672        endif
 673       endif
 674       if(me.eq.king.or..not.out1file)
 675      & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
 676       return
 677       end
 678 c------------------------------------------------------------------------------
 679       subroutine molread
 680 C
 681 C Read molecular data.
 682 C
 683       implicit none
 684       include 'DIMENSIONS'
 685 #ifdef MPI
 686       include 'mpif.h'
 687       integer error_msg,ierror,ierr,ierrcode
 688 #endif
 689       include 'COMMON.IOUNITS'
 690       include 'COMMON.GEO'
 691       include 'COMMON.VAR'
 692       include 'COMMON.INTERACT'
 693       include 'COMMON.LOCAL'
 694       include 'COMMON.NAMES'
 695       include 'COMMON.CHAIN'
 696       include 'COMMON.FFIELD'
 697       include 'COMMON.SBRIDGE'
 698       include 'COMMON.HEADER'
 699       include 'COMMON.CONTROL'
 700       include 'COMMON.SAXS'
 701       include 'COMMON.DBASE'
 702       include 'COMMON.THREAD'
 703       include 'COMMON.CONTACTS'
 704       include 'COMMON.TORCNSTR'
 705       include 'COMMON.TIME1'
 706       include 'COMMON.BOUNDS'
 707       include 'COMMON.MD'
 708       include 'COMMON.SETUP'
 709       include 'COMMON.SHIELD'
 710       character*4 sequence(maxres)
 711       integer rescode
 712       double precision x(maxvar)
 713       character*256 pdbfile
 714       character*400 weightcard
 715       character*80 weightcard_t,ucase
 716       integer itype_pdb(maxres)
 717       common /pizda/ itype_pdb
 718       logical seq_comp,fail
 719       double precision energia(0:n_ene)
 720       double precision secprob(3,maxdih_constr)
 721       double precision co
 722       double precision phihel,phibet,sigmahel,sigmabet
 723       integer iti,nsi,maxsi
 724       integer ilen
 725       external ilen
 726       integer iperm,tperm
 727       integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2
 728       double precision sumv
 729 C
 730 C Read PDB structure if applicable
 731 C
 732       if (indpdb.gt.0 .or. pdbref) then
 733         read(inp,'(a)') pdbfile
 734         if(me.eq.king.or..not.out1file)
 735      &   write (iout,'(2a)') 'PDB data will be read from file ',
 736      &   pdbfile(:ilen(pdbfile))
 737         open(ipdbin,file=pdbfile,status='old',err=33)
 738         goto 34
 739   33    write (iout,'(a)') 'Error opening PDB file.'
 740         stop
 741   34    continue
 742         call readpdb
 743         do i=1,2*nres
 744           do j=1,3
 745             crefjlee(j,i)=c(j,i)
 746           enddo
 747         enddo
 748 #ifdef DEBUG
 749         do i=1,nres
 750           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
 751      &      (crefjlee(j,i+nres),j=1,3)
 752         enddo
 753 #endif
 754         if(me.eq.king.or..not.out1file)
 755      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 756      &   ' nstart_sup=',nstart_sup
 757         do i=1,nres
 758           itype_pdb(i)=itype(i)
 759         enddo
 760         close (ipdbin)
 761         nnt=nstart_sup
 762         nct=nstart_sup+nsup-1
 763         call contact(.false.,ncont_ref,icont_ref,co)
 764
 765         if (sideadd) then
 766          if(me.eq.king.or..not.out1file)
 767      &    write(iout,*)'Adding sidechains'
 768          maxsi=1000
 769          do i=2,nres-1
 770           iti=itype(i)
 771           if (iti.ne.10 .and. itype(i).ne.ntyp1) then
 772             nsi=0
 773             fail=.true.
 774             do while (fail.and.nsi.le.maxsi)
 775               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 776               nsi=nsi+1
 777             enddo
 778             if(fail) write(iout,*)'Adding sidechain failed for res ',
 779      &              i,' after ',nsi,' trials'
 780           endif
 781          enddo
 782         endif
 783       endif
 784       if (indpdb.eq.0) then
 785 C Read sequence if not taken from the pdb file.
 786         read (inp,*) nres
 787 c        print *,'nres=',nres
 788         if (iscode.gt.0) then
 789           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 790         else
 791           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 792         endif
 793 C Convert sequence to numeric code
 794         do i=1,nres
 795           itype(i)=rescode(i,sequence(i),iscode)
 796         enddo
 797       endif
 798 c      print *,nres
 799 c      print '(20i4)',(itype(i),i=1,nres)
 800       do i=1,nres
 801 #ifdef PROCOR
 802         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 803 #else
 804         if (itype(i).eq.ntyp1) then
 805 #endif
 806           itel(i)=0
 807 #ifdef PROCOR
 808         else if (iabs(itype(i+1)).ne.20) then
 809 #else
 810         else if (iabs(itype(i)).ne.20) then
 811 #endif
 812           itel(i)=1
 813         else
 814           itel(i)=2
 815         endif
 816       enddo
 817       if(me.eq.king.or..not.out1file)then
 818        write (iout,*) "ITEL"
 819        do i=1,nres-1
 820          write (iout,*) i,itype(i),itel(i)
 821        enddo
 822        print *,'Call Read_Bridge.'
 823       endif
 824       nnt=1
 825       nct=nres
 826 cd      print *,'NNT=',NNT,' NCT=',NCT
 827       call seq2chains(nres,itype,nchain,chain_length,chain_border,
 828      &  ireschain)
 829       chain_border1(1,1)=1
 830       chain_border1(2,1)=chain_border(2,1)+1
 831       do i=2,nchain-1
 832         chain_border1(1,i)=chain_border(1,i)-1
 833         chain_border1(2,i)=chain_border(2,i)+1
 834       enddo
 835       chain_border1(1,nchain)=chain_border(1,nchain)-1
 836       chain_border1(2,nchain)=nres
 837       write(iout,*) "nres",nres," nchain",nchain
 838       do i=1,nchain
 839         write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
 840      &    chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
 841       enddo
 842       call chain_symmetry(nchain,nres,itype,chain_border,
 843      &    chain_length,npermchain,tabpermchain)
 844 c      do i=1,nres
 845 c        write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
 846 c     &    ii=1,npermchain)
 847 c      enddo
 848       write(iout,*) "residue permutations"
 849       do i=1,nres
 850         write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
 851       enddo
 852       call flush(iout)
 853       if (itype(1).eq.ntyp1) nnt=2
 854       if (itype(nres).eq.ntyp1) nct=nct-1
 855       write (iout,*) "nnt",nnt," nct",nct
 856       call flush(iout)
 857 #ifdef DFA
 858       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
 859      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
 860        call init_dfa_vars
 861        print*, 'init_dfa_vars finished!'
 862        call read_dfa_info
 863        print*, 'read_dfa_info finished!'
 864       endif
 865 #endif
 866       if (pdbref) then
 867         if(me.eq.king.or..not.out1file)
 868      &   write (iout,'(a,i3)') 'nsup=',nsup
 869         nstart_seq=nnt
 870         if (nsup.le.(nct-nnt+1)) then
 871           do i=0,nct-nnt+1-nsup
 872             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 873               nstart_seq=nnt+i
 874               goto 111
 875             endif
 876           enddo
 877           write (iout,'(a)')
 878      &            'Error - sequences to be superposed do not match.'
 879           stop
 880         else
 881           do i=0,nsup-(nct-nnt+1)
 882             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 883      &      then
 884               nstart_sup=nstart_sup+i
 885               nsup=nct-nnt+1
 886               goto 111
 887             endif
 888           enddo
 889           write (iout,'(a)')
 890      &            'Error - sequences to be superposed do not match.'
 891         endif
 892   111   continue
 893         if (nsup.eq.0) nsup=nct-nnt
 894         if (nstart_sup.eq.0) nstart_sup=nnt
 895         if (nstart_seq.eq.0) nstart_seq=nnt
 896         if(me.eq.king.or..not.out1file)
 897      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 898      &                 ' nstart_seq=',nstart_seq
 899       endif
 900       call read_bridge
 901 C 8/13/98 Set limits to generating the dihedral angles
 902       do i=1,nres
 903         phibound(1,i)=-pi
 904         phibound(2,i)=pi
 905       enddo
 906       read (inp,*) ndih_constr
 907       write (iout,*) "ndih_constr",ndih_constr
 908       if (ndih_constr.gt.0) then
 909         raw_psipred=.false.
 910 C        read (inp,*) ftors
 911         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
 912      &  i=1,ndih_constr)
 913 #ifdef MPI
 914         if(me.eq.king.or..not.out1file)then
 915 #endif
 916          write (iout,*)
 917      &   'There are',ndih_constr,' restraints on gamma angles.'
 918          do i=1,ndih_constr
 919           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
 920      &    ftors(i)
 921          enddo
 922 #ifdef MPI
 923         endif
 924 #endif
 925         do i=1,ndih_constr
 926           phi0(i)=deg2rad*phi0(i)
 927           drange(i)=deg2rad*drange(i)
 928         enddo
 929 C        if(me.eq.king.or..not.out1file)
 930 C     &   write (iout,*) 'FTORS',ftors
 931         do i=1,ndih_constr
 932           ii = idih_constr(i)
 933           phibound(1,ii) = phi0(i)-drange(i)
 934           phibound(2,ii) = phi0(i)+drange(i)
 935         enddo
 936 #ifdef MPI
 937         if (me.eq.king .or. .not.out1file) then
 938 #endif
 939          write (iout,'(a)') 'Boundaries in gamma angle sampling:'
 940          do i=1,nres
 941            write (iout,'(a3,i5,2f10.1)')
 942      &    restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 943          enddo
 944 #ifdef MP
 945         endif
 946 #endif
 947       else if (ndih_constr.lt.0) then
 948         raw_psipred=.true.
 949         call card_concat(weightcard)
 950         call reada(weightcard,"PHIHEL",phihel,50.0D0)
 951         call reada(weightcard,"PHIBET",phibet,180.0D0)
 952         call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
 953         call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
 954         call reada(weightcard,"WDIHC",wdihc,0.591D0)
 955         write (iout,*) "Weight of dihedral angle restraints",wdihc
 956         read(inp,'(9x,3f7.3)')
 957      &     (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
 958         write (iout,*) "The secprob array"
 959         do i=nnt,nct
 960           write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
 961         enddo
 962         ndih_constr=0
 963         do i=nnt+3,nct
 964           if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
 965      &    .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
 966             ndih_constr=ndih_constr+1
 967             idih_constr(ndih_constr)=i
 968             sumv=0.0d0
 969             do j=1,3
 970               vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
 971               sumv=sumv+vpsipred(j,ndih_constr)
 972             enddo
 973             do j=1,3
 974               vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
 975             enddo
 976             phibound(1,ndih_constr)=phihel*deg2rad
 977             phibound(2,ndih_constr)=phibet*deg2rad
 978             sdihed(1,ndih_constr)=sigmahel*deg2rad
 979             sdihed(2,ndih_constr)=sigmabet*deg2rad
 980           endif
 981         enddo
 982 #ifdef MPI
 983         if(me.eq.king.or..not.out1file)then
 984 #endif
 985          write (iout,*)
 986      &   'There are',ndih_constr,
 987      &   ' bimodal restraints on gamma angles.'
 988          do i=1,ndih_constr
 989            write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
 990      &      restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
 991      &      restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
 992      &      phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
 993      &      phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
 994      &      (vpsipred(j,i),j=1,3)
 995          enddo
 996 #ifdef MPI
 997         endif
 998 #endif
 999 c Raw psipred input
1000       endif
1001 C first setting the theta boundaries to 0 to pi
1002 C this mean that there is no energy penalty for any angle occuring this can be applied
1003 C for generate random conformation but is not implemented in this way
1004 C      do i=1,nres
1005 C        thetabound(1,i)=0
1006 C        thetabound(2,i)=pi
1007 C      enddo
1008 C begin reading theta constrains this is quartic constrains allowing to
1009 C have smooth second derivative
1010       if (with_theta_constr) then
1011 C with_theta_constr is keyword allowing for occurance of theta constrains
1012       read (inp,*) ntheta_constr
1013 C ntheta_constr is the number of theta constrains
1014       if (ntheta_constr.gt.0) then
1015 C        read (inp,*) ftors
1016         read (inp,*) (itheta_constr(i),theta_constr0(i),
1017      &  theta_drange(i),for_thet_constr(i),
1018      &  i=1,ntheta_constr)
1019 C the above code reads from 1 to ntheta_constr
1020 C itheta_constr(i) residue i for which is theta_constr
1021 C theta_constr0 the global minimum value
1022 C theta_drange is range for which there is no energy penalty
1023 C for_thet_constr is the force constant for quartic energy penalty
1024 C E=k*x**4
1025         if(me.eq.king.or..not.out1file)then
1026          write (iout,*)
1027      &   'There are',ntheta_constr,' constraints on phi angles.'
1028          do i=1,ntheta_constr
1029           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1030      &    theta_drange(i),
1031      &    for_thet_constr(i)
1032          enddo
1033         endif
1034         do i=1,ntheta_constr
1035           theta_constr0(i)=deg2rad*theta_constr0(i)
1036           theta_drange(i)=deg2rad*theta_drange(i)
1037         enddo
1038 C        if(me.eq.king.or..not.out1file)
1039 C     &   write (iout,*) 'FTORS',ftors
1040 C        do i=1,ntheta_constr
1041 C          ii = itheta_constr(i)
1042 C          thetabound(1,ii) = phi0(i)-drange(i)
1043 C          thetabound(2,ii) = phi0(i)+drange(i)
1044 C        enddo
1045       endif ! ntheta_constr.gt.0
1046       endif! with_theta_constr
1047 C
1048 C      with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1049 C      write (iout,*) "with_dihed_constr ",with_dihed_constr
1050 c--- Zscore rms -------
1051       if (nz_start.eq.0) nz_start=nnt
1052       if (nz_end.eq.0 .and. nsup.gt.0) then
1053         nz_end=nnt+nsup-1
1054       else if (nz_end.eq.0) then
1055         nz_end=nct
1056       endif
1057       if(me.eq.king.or..not.out1file)then
1058        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1059        write (iout,*) 'IZ_SC=',iz_sc
1060       endif
1061 c----------------------
1062       call init_int_table
1063       if (refstr) then
1064         if (.not.pdbref) then
1065           call read_angles(inp,*38)
1066           call bond_regular
1067           goto 39
1068    38     write (iout,'(a)') 'Error reading reference structure.'
1069 #ifdef MPI
1070           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1071           stop 'Error reading reference structure'
1072 #endif
1073    39     call chainbuild_extconf
1074           call setup_var
1075 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1076           nstart_sup=nnt
1077           nstart_seq=nnt
1078           nsup=nct-nnt+1
1079           do i=1,2*nres
1080             do j=1,3
1081               cref(j,i)=c(j,i)
1082             enddo
1083           enddo
1084           call contact(.true.,ncont_ref,icont_ref,co)
1085         endif
1086       endif
1087 c        print *, "A TU"
1088       write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1089       call flush(iout)
1090       if (constr_dist.gt.0) call read_dist_constr
1091       write (iout,*) "After read_dist_constr nhpb",nhpb
1092       if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1093       call hpb_partition
1094       call NMRpeak_partition
1095       if(me.eq.king.or..not.out1file)write (iout,*) 'Contact order:',co
1096       if (pdbref) then
1097         if(me.eq.king.or..not.out1file)
1098      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1099         do i=1,ncont_ref
1100           do j=1,2
1101             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1102           enddo
1103           if(me.eq.king.or..not.out1file)
1104      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1105      &     icont_ref(1,i),' ',
1106      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1107         enddo
1108       endif
1109       write (iout,*) "calling read_saxs_consrtr",nsaxs
1110       if (nsaxs.gt.0) call read_saxs_constr
1111
1112       if (constr_homology.gt.0) then
1113         call read_constr_homology
1114         if (indpdb.gt.0 .or. pdbref) then
1115           do i=1,2*nres
1116             do j=1,3
1117               c(j,i)=crefjlee(j,i)
1118               cref(j,i)=crefjlee(j,i)
1119             enddo
1120           enddo
1121         endif
1122 #ifdef DEBUG
1123         write (iout,*) "sc_loc_geom: Array C"
1124         do i=1,nres
1125           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1126      &      (c(j,i+nres),j=1,3)
1127         enddo
1128         write (iout,*) "Array Cref"
1129         do i=1,nres
1130           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1131      &      (cref(j,i+nres),j=1,3)
1132         enddo
1133 #endif
1134        call int_from_cart1(.false.)
1135        call sc_loc_geom(.false.)
1136        do i=1,nres
1137          thetaref(i)=theta(i)
1138          phiref(i)=phi(i)
1139        enddo
1140        do i=1,nres-1
1141          do j=1,3
1142            dc(j,i)=c(j,i+1)-c(j,i)
1143            dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1144          enddo
1145        enddo
1146        do i=2,nres-1
1147          do j=1,3
1148            dc(j,i+nres)=c(j,i+nres)-c(j,i)
1149            dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1150          enddo
1151        enddo
1152       else
1153         homol_nset=0
1154       endif
1155
1156
1157 C      endif
1158       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1159      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1160      &    modecalc.ne.10) then
1161 C If input structure hasn't been supplied from the PDB file read or generate
1162 C initial geometry.
1163         if (iranconf.eq.0 .and. .not. extconf) then
1164           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1165      &     write (iout,'(a)') 'Initial geometry will be read in.'
1166           if (read_cart) then
1167             read(inp,'(8f10.5)',end=36,err=36)
1168      &       ((c(l,k),l=1,3),k=1,nres),
1169      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1170             write (iout,*) "Exit READ_CART"
1171 c            write (iout,'(8f10.5)')
1172 c     &       ((c(l,k),l=1,3),k=1,nres),
1173 c     &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1174             call cartprint
1175             do i=1,nres-1
1176               do j=1,3
1177                 dc(j,i)=c(j,i+1)-c(j,i)
1178 c                dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1179               enddo
1180             enddo
1181             do i=nnt,nct
1182               if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1183                 do j=1,3
1184                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1185 c                  dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1186                 enddo
1187               else
1188                 do j=1,3
1189                   dc(j,i+nres)=0.0d0
1190 c                  dc_norm(j,i+nres)=0.0d0
1191                 enddo
1192               endif
1193             enddo
1194             call int_from_cart1(.true.)
1195             write (iout,*) "Finish INT_TO_CART"
1196 c      write (iout,*) "DC"
1197 c      do i=1,nres
1198 c        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1199 c     &    (dc(j,i+nres),j=1,3)
1200 c      enddo
1201 c            call cartprint
1202 c            call setup_var
1203 c            return
1204           else
1205             call read_angles(inp,*36)
1206             call bond_regular
1207             call chainbuild_extconf
1208           endif
1209           goto 37
1210    36     write (iout,'(a)') 'Error reading angle file.'
1211 #ifdef MPI
1212           call mpi_finalize( MPI_COMM_WORLD,IERR )
1213 #endif
1214           stop 'Error reading angle file.'
1215    37     continue
1216         else if (extconf) then
1217           if (me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1218      &      write (iout,'(a)') 'Extended chain initial geometry.'
1219           do i=3,nres
1220             theta(i)=90d0*deg2rad
1221           enddo
1222           do i=4,nres
1223             phi(i)=180d0*deg2rad
1224           enddo
1225           do i=2,nres-1
1226             alph(i)=110d0*deg2rad
1227           enddo
1228           do i=2,nres-1
1229             omeg(i)=-120d0*deg2rad
1230             if (itype(i).le.0) omeg(i)=-omeg(i)
1231           enddo
1232           call bond_regular
1233           call chainbuild_extconf
1234         else
1235           if(me.eq.king.or..not.out1file)
1236      &     write (iout,'(a)') 'Random-generated initial geometry.'
1237           call bond_regular
1238 #ifdef MPI
1239           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1240      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1241 #endif
1242             do itrial=1,100
1243               itmp=1
1244               call gen_rand_conf(itmp,*30)
1245               goto 40
1246    30         write (iout,*) 'Failed to generate random conformation',
1247      &          ', itrial=',itrial
1248               write (*,*) 'Processor:',me,
1249      &          ' Failed to generate random conformation',
1250      &          ' itrial=',itrial
1251               call intout
1252
1253 #ifdef AIX
1254               call flush_(iout)
1255 #else
1256               call flush(iout)
1257 #endif
1258             enddo
1259             write (iout,'(a,i3,a)') 'Processor:',me,
1260      &        ' error in generating random conformation.'
1261             write (*,'(a,i3,a)') 'Processor:',me,
1262      &        ' error in generating random conformation.'
1263             call flush(iout)
1264 #ifdef MPI
1265             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1266    40       continue
1267           endif
1268 #else
1269           write (*,'(a)')
1270      &      ' error in generating random conformation.'
1271           stop
1272    40     continue
1273 #endif
1274         endif
1275       elseif (modecalc.eq.4) then
1276         read (inp,'(a)') intinname
1277         open (intin,file=intinname,status='old',err=333)
1278         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1279      &  write (iout,'(a)') 'intinname',intinname
1280         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1281         goto 334
1282   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1283 #ifdef MPI
1284         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1285 #endif
1286         stop 'Error opening angle file.'
1287   334   continue
1288
1289       endif
1290 C Generate distance constraints, if the PDB structure is to be regularized.
1291       if (nthread.gt.0) then
1292         call read_threadbase
1293       endif
1294       call setup_var
1295       if (me.eq.king .or. .not. out1file)
1296      & call intout
1297       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1298         write (iout,'(/a,i3,a)')
1299      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1300         write (iout,'(20i4)') (iss(i),i=1,ns)
1301        if (dyn_ss) then
1302           write(iout,*)"Running with dynamic disulfide-bond formation"
1303        else
1304         write (iout,'(/a/)') 'Pre-formed links are:'
1305         do i=1,nss
1306           i1=ihpb(i)-nres
1307           i2=jhpb(i)-nres
1308           it1=itype(i1)
1309           it2=itype(i2)
1310           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1311      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1312      &    ebr,forcon(i)
1313         enddo
1314         write (iout,'(a)')
1315        endif
1316       endif
1317       if (ns.gt.0.and.dyn_ss) then
1318           do i=nss+1,nhpb
1319             ihpb(i-nss)=ihpb(i)
1320             jhpb(i-nss)=jhpb(i)
1321             forcon(i-nss)=forcon(i)
1322             dhpb(i-nss)=dhpb(i)
1323           enddo
1324           nhpb=nhpb-nss
1325           nss=0
1326           call hpb_partition
1327           do i=1,ns
1328             dyn_ss_mask(iss(i))=.true.
1329           enddo
1330       endif
1331 c --- test
1332 c      call cartprint
1333 c      write (iout,*) "DC"
1334 c      do i=1,nres
1335 c        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1336 c     &    (dc(j,i+nres),j=1,3)
1337 c      enddo
1338       if (i2ndstr.gt.0) call secstrp2dihc
1339 c      call geom_to_var(nvar,x)
1340 c      call etotal(energia(0))
1341 c      call enerprint(energia(0))
1342 c      call briefout(0,etot)
1343 c      stop
1344 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1345 cd    write (iout,'(a)') 'Variable list:'
1346 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1347 #ifdef MPI
1348       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1349      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1350      &  'Processor',myrank,': end reading molecular data.'
1351 #endif
1352 c      print *,"A TU?"
1353       return
1354       end
1355 c--------------------------------------------------------------------------
1356       logical function seq_comp(itypea,itypeb,length)
1357       implicit none
1358       integer length,itypea(length),itypeb(length)
1359       integer i
1360       do i=1,length
1361         if (itypea(i).ne.itypeb(i)) then
1362           seq_comp=.false.
1363           return
1364         endif
1365       enddo
1366       seq_comp=.true.
1367       return
1368       end
1369 c-----------------------------------------------------------------------------
1370       subroutine read_bridge
1371 C Read information about disulfide bridges.
1372       implicit none
1373       include 'DIMENSIONS'
1374 #ifdef MPI
1375       include 'mpif.h'
1376       integer ierror
1377 #endif
1378       include 'COMMON.IOUNITS'
1379       include 'COMMON.GEO'
1380       include 'COMMON.VAR'
1381       include 'COMMON.INTERACT'
1382       include 'COMMON.LOCAL'
1383       include 'COMMON.NAMES'
1384       include 'COMMON.CHAIN'
1385       include 'COMMON.FFIELD'
1386       include 'COMMON.SBRIDGE'
1387       include 'COMMON.HEADER'
1388       include 'COMMON.CONTROL'
1389       include 'COMMON.DBASE'
1390       include 'COMMON.THREAD'
1391       include 'COMMON.TIME1'
1392       include 'COMMON.SETUP'
1393       integer i,j
1394 C Read bridging residues.
1395       read (inp,*) ns,(iss(i),i=1,ns)
1396       print *,'ns=',ns
1397       if(me.eq.king.or..not.out1file)
1398      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1399 C Check whether the specified bridging residues are cystines.
1400       do i=1,ns
1401         if (itype(iss(i)).ne.1) then
1402           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1403      &   'Do you REALLY think that the residue ',
1404      &    restyp(itype(iss(i))),i,
1405      &   ' can form a disulfide bridge?!!!'
1406           write (*,'(2a,i3,a)')
1407      &   'Do you REALLY think that the residue ',
1408      &    restyp(itype(iss(i))),i,
1409      &   ' can form a disulfide bridge?!!!'
1410 #ifdef MPI
1411          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1412          stop
1413 #endif
1414         endif
1415       enddo
1416 C Read preformed bridges.
1417       if (ns.gt.0) then
1418       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1419       if(fg_rank.eq.0)
1420      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1421       if (nss.gt.0) then
1422         nhpb=nss
1423 C Check if the residues involved in bridges are in the specified list of
1424 C bridging residues.
1425         do i=1,nss
1426           do j=1,i-1
1427             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1428      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1429               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1430      &      ' contains residues present in other pairs.'
1431               write (*,'(a,i3,a)') 'Disulfide pair',i,
1432      &      ' contains residues present in other pairs.'
1433 #ifdef MPI
1434               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1435               stop
1436 #endif
1437             endif
1438           enddo
1439           do j=1,ns
1440             if (ihpb(i).eq.iss(j)) goto 10
1441           enddo
1442           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1443    10     continue
1444           do j=1,ns
1445             if (jhpb(i).eq.iss(j)) goto 20
1446           enddo
1447           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1448    20     continue
1449 c          dhpb(i)=dbr
1450 c          forcon(i)=fbr
1451         enddo
1452         do i=1,nss
1453           ihpb(i)=ihpb(i)+nres
1454           jhpb(i)=jhpb(i)+nres
1455         enddo
1456       endif
1457       endif
1458       return
1459       end
1460 c----------------------------------------------------------------------------
1461       subroutine read_x(kanal,*)
1462       implicit none
1463       include 'DIMENSIONS'
1464       include 'COMMON.GEO'
1465       include 'COMMON.VAR'
1466       include 'COMMON.CHAIN'
1467       include 'COMMON.IOUNITS'
1468       include 'COMMON.CONTROL'
1469       include 'COMMON.LOCAL'
1470       include 'COMMON.INTERACT'
1471       integer i,j,k,l,kanal
1472 c Read coordinates from input
1473 c
1474       read(kanal,'(8f10.5)',end=10,err=10)
1475      &  ((c(l,k),l=1,3),k=1,nres),
1476      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1477       do j=1,3
1478         c(j,nres+1)=c(j,1)
1479         c(j,2*nres)=c(j,nres)
1480       enddo
1481       call int_from_cart1(.false.)
1482       do i=1,nres-1
1483         do j=1,3
1484           dc(j,i)=c(j,i+1)-c(j,i)
1485           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1486         enddo
1487       enddo
1488       do i=nnt,nct
1489         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1490           do j=1,3
1491             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1492             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1493           enddo
1494         endif
1495       enddo
1496
1497       return
1498    10 return1
1499       end
1500 c----------------------------------------------------------------------------
1501       subroutine read_threadbase
1502       implicit none
1503       include 'DIMENSIONS'
1504       include 'COMMON.IOUNITS'
1505       include 'COMMON.GEO'
1506       include 'COMMON.VAR'
1507       include 'COMMON.INTERACT'
1508       include 'COMMON.LOCAL'
1509       include 'COMMON.NAMES'
1510       include 'COMMON.CHAIN'
1511       include 'COMMON.FFIELD'
1512       include 'COMMON.SBRIDGE'
1513       include 'COMMON.HEADER'
1514       include 'COMMON.CONTROL'
1515       include 'COMMON.DBASE'
1516       include 'COMMON.THREAD'
1517       include 'COMMON.TIME1'
1518       integer i,j,k
1519       double precision dist
1520 C Read pattern database for threading.
1521       read (icbase,*) nseq
1522       do i=1,nseq
1523         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1524      &   nres_base(2,i),nres_base(3,i)
1525         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1526      &   nres_base(1,i))
1527 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1528 c    &   nres_base(2,i),nres_base(3,i)
1529 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1530 c    &   nres_base(1,i))
1531       enddo
1532       close (icbase)
1533       if (weidis.eq.0.0D0) weidis=0.1D0
1534       do i=nnt,nct
1535         do j=i+2,nct
1536           nhpb=nhpb+1
1537           ihpb(nhpb)=i
1538           jhpb(nhpb)=j
1539           forcon(nhpb)=weidis
1540         enddo
1541       enddo
1542       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1543       write (iout,'(a,i5)') 'nexcl: ',nexcl
1544       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1545       return
1546       end
1547 c------------------------------------------------------------------------------
1548       subroutine setup_var
1549       implicit none
1550       include 'DIMENSIONS'
1551       include 'COMMON.IOUNITS'
1552       include 'COMMON.GEO'
1553       include 'COMMON.VAR'
1554       include 'COMMON.INTERACT'
1555       include 'COMMON.LOCAL'
1556       include 'COMMON.NAMES'
1557       include 'COMMON.CHAIN'
1558       include 'COMMON.FFIELD'
1559       include 'COMMON.SBRIDGE'
1560       include 'COMMON.HEADER'
1561       include 'COMMON.CONTROL'
1562       include 'COMMON.DBASE'
1563       include 'COMMON.THREAD'
1564       include 'COMMON.TIME1'
1565       integer i
1566 C Set up variable list.
1567       ntheta=nres-2
1568       nphi=nres-3
1569       nvar=ntheta+nphi
1570       nside=0
1571       do i=2,nres-1
1572         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1573           nside=nside+1
1574           ialph(i,1)=nvar+nside
1575           ialph(nside,2)=i
1576         endif
1577       enddo
1578       if (indphi.gt.0) then
1579         nvar=nphi
1580       else if (indback.gt.0) then
1581         nvar=nphi+ntheta
1582       else
1583         nvar=nvar+2*nside
1584       endif
1585       return
1586       end
1587 c----------------------------------------------------------------------------
1588       subroutine gen_dist_constr
1589 C Generate CA distance constraints.
1590       implicit none
1591       include 'DIMENSIONS'
1592       include 'COMMON.IOUNITS'
1593       include 'COMMON.GEO'
1594       include 'COMMON.VAR'
1595       include 'COMMON.INTERACT'
1596       include 'COMMON.LOCAL'
1597       include 'COMMON.NAMES'
1598       include 'COMMON.CHAIN'
1599       include 'COMMON.FFIELD'
1600       include 'COMMON.SBRIDGE'
1601       include 'COMMON.HEADER'
1602       include 'COMMON.CONTROL'
1603       include 'COMMON.DBASE'
1604       include 'COMMON.THREAD'
1605       include 'COMMON.TIME1'
1606       integer i,j,itype_pdb(maxres)
1607       common /pizda/ itype_pdb
1608       double precision dist
1609       character*2 iden
1610 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1611 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1612 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1613 cd     & ' nsup',nsup
1614       do i=nstart_sup,nstart_sup+nsup-1
1615 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1616 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1617         do j=i+2,nstart_sup+nsup-1
1618           nhpb=nhpb+1
1619           ihpb(nhpb)=i+nstart_seq-nstart_sup
1620           jhpb(nhpb)=j+nstart_seq-nstart_sup
1621           forcon(nhpb)=weidis
1622           dhpb(nhpb)=dist(i,j)
1623         enddo
1624       enddo
1625 cd      write (iout,'(a)') 'Distance constraints:'
1626 cd      do i=nss+1,nhpb
1627 cd        ii=ihpb(i)
1628 cd        jj=jhpb(i)
1629 cd        iden='CA'
1630 cd        if (ii.gt.nres) then
1631 cd          iden='SC'
1632 cd          ii=ii-nres
1633 cd          jj=jj-nres
1634 cd        endif
1635 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1636 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1637 cd     &  dhpb(i),forcon(i)
1638 cd      enddo
1639       return
1640       end
1641 c----------------------------------------------------------------------------
1642       subroutine map_read
1643       implicit none
1644       include 'DIMENSIONS'
1645       include 'COMMON.MAP'
1646       include 'COMMON.IOUNITS'
1647       integer imap
1648       character*3 angid(4) /'THE','PHI','ALP','OME'/
1649       character*80 mapcard,ucase
1650       do imap=1,nmap
1651         read (inp,'(a)') mapcard
1652         mapcard=ucase(mapcard)
1653         if (index(mapcard,'PHI').gt.0) then
1654           kang(imap)=1
1655         else if (index(mapcard,'THE').gt.0) then
1656           kang(imap)=2
1657         else if (index(mapcard,'ALP').gt.0) then
1658           kang(imap)=3
1659         else if (index(mapcard,'OME').gt.0) then
1660           kang(imap)=4
1661         else
1662           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1663           stop 'Error - illegal variable spec in MAP card.'
1664         endif
1665         call readi (mapcard,'RES1',res1(imap),0)
1666         call readi (mapcard,'RES2',res2(imap),0)
1667         if (res1(imap).eq.0) then
1668           res1(imap)=res2(imap)
1669         else if (res2(imap).eq.0) then
1670           res2(imap)=res1(imap)
1671         endif
1672         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1673           write (iout,'(a)')
1674      &    'Error - illegal definition of variable group in MAP.'
1675           stop 'Error - illegal definition of variable group in MAP.'
1676         endif
1677         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1678         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1679         call readi(mapcard,'NSTEP',nstep(imap),0)
1680         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1681           write (iout,'(a)')
1682      &     'Illegal boundary and/or step size specification in MAP.'
1683           stop 'Illegal boundary and/or step size specification in MAP.'
1684         endif
1685       enddo ! imap
1686       return
1687       end
1688 c----------------------------------------------------------------------------
1689       subroutine csaread
1690       implicit none
1691       include 'DIMENSIONS'
1692       include 'COMMON.IOUNITS'
1693       include 'COMMON.GEO'
1694       include 'COMMON.CSA'
1695       include 'COMMON.BANK'
1696       include 'COMMON.CONTROL'
1697       character*80 ucase
1698       character*620 mcmcard
1699       call card_concat(mcmcard)
1700
1701       call readi(mcmcard,'NCONF',nconf,50)
1702       call readi(mcmcard,'NADD',nadd,0)
1703       call readi(mcmcard,'JSTART',jstart,1)
1704       call readi(mcmcard,'JEND',jend,1)
1705       call readi(mcmcard,'NSTMAX',nstmax,500000)
1706       call readi(mcmcard,'N0',n0,1)
1707       call readi(mcmcard,'N1',n1,6)
1708       call readi(mcmcard,'N2',n2,4)
1709       call readi(mcmcard,'N3',n3,0)
1710       call readi(mcmcard,'N4',n4,0)
1711       call readi(mcmcard,'N5',n5,0)
1712       call readi(mcmcard,'N6',n6,10)
1713       call readi(mcmcard,'N7',n7,0)
1714       call readi(mcmcard,'N8',n8,0)
1715       call readi(mcmcard,'N9',n9,0)
1716       call readi(mcmcard,'N14',n14,0)
1717       call readi(mcmcard,'N15',n15,0)
1718       call readi(mcmcard,'N16',n16,0)
1719       call readi(mcmcard,'N17',n17,0)
1720       call readi(mcmcard,'N18',n18,0)
1721
1722       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1723
1724       call readi(mcmcard,'NDIFF',ndiff,2)
1725       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1726       call readi(mcmcard,'IS1',is1,1)
1727       call readi(mcmcard,'IS2',is2,8)
1728       call readi(mcmcard,'NRAN0',nran0,4)
1729       call readi(mcmcard,'NRAN1',nran1,2)
1730       call readi(mcmcard,'IRR',irr,1)
1731       call readi(mcmcard,'NSEED',nseed,20)
1732       call readi(mcmcard,'NTOTAL',ntotal,10000)
1733       call reada(mcmcard,'CUT1',cut1,2.0d0)
1734       call reada(mcmcard,'CUT2',cut2,5.0d0)
1735       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1736       call readi(mcmcard,'ICMAX',icmax,3)
1737       call readi(mcmcard,'IRESTART',irestart,0)
1738 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1739       ntbankm=0
1740 c!bankt
1741       call reada(mcmcard,'DELE',dele,20.0d0)
1742       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1743       call readi(mcmcard,'IREF',iref,0)
1744       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1745       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1746       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1747       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1748       write (iout,*) "NCONF_IN",nconf_in
1749       return
1750       end
1751 c----------------------------------------------------------------------------
1752       subroutine mcmread
1753       implicit none
1754       include 'DIMENSIONS'
1755       include 'COMMON.MCM'
1756       include 'COMMON.MCE'
1757       include 'COMMON.IOUNITS'
1758       character*80 ucase
1759       character*320 mcmcard
1760       integer i
1761       call card_concat(mcmcard)
1762       call readi(mcmcard,'MAXACC',maxacc,100)
1763       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1764       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1765       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1766       call readi(mcmcard,'MAXREPM',maxrepm,200)
1767       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1768       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1769       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1770       call reada(mcmcard,'E_UP',e_up,5.0D0)
1771       call reada(mcmcard,'DELTE',delte,0.1D0)
1772       call readi(mcmcard,'NSWEEP',nsweep,5)
1773       call readi(mcmcard,'NSTEPH',nsteph,0)
1774       call readi(mcmcard,'NSTEPC',nstepc,0)
1775       call reada(mcmcard,'TMIN',tmin,298.0D0)
1776       call reada(mcmcard,'TMAX',tmax,298.0D0)
1777       call readi(mcmcard,'NWINDOW',nwindow,0)
1778       call readi(mcmcard,'PRINT_MC',print_mc,0)
1779       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1780       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1781       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1782       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1783       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1784       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1785       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1786       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1787       if (nwindow.gt.0) then
1788         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1789         do i=1,nwindow
1790           winlen(i)=winend(i)-winstart(i)+1
1791         enddo
1792       endif
1793       if (tmax.lt.tmin) tmax=tmin
1794       if (tmax.eq.tmin) then
1795         nstepc=0
1796         nsteph=0
1797       endif
1798       if (nstepc.gt.0 .and. nsteph.gt.0) then
1799         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1800         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1801       endif
1802 C Probabilities of different move types
1803       sumpro_type(0)=0.0D0
1804       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1805       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1806       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1807       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1808       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1809       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1810       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1811       do i=1,MaxMoveType
1812         print *,'i',i,' sumprotype',sumpro_type(i)
1813         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1814         print *,'i',i,' sumprotype',sumpro_type(i)
1815       enddo
1816       return
1817       end
1818 c----------------------------------------------------------------------------
1819       subroutine read_minim
1820       implicit none
1821       include 'DIMENSIONS'
1822       include 'COMMON.MINIM'
1823       include 'COMMON.IOUNITS'
1824       character*80 ucase
1825       character*320 minimcard
1826       call card_concat(minimcard)
1827       call readi(minimcard,'MAXMIN',maxmin,2000)
1828       call readi(minimcard,'MAXFUN',maxfun,5000)
1829       call readi(minimcard,'MINMIN',minmin,maxmin)
1830       call readi(minimcard,'MINFUN',minfun,maxmin)
1831       call reada(minimcard,'TOLF',tolf,1.0D-2)
1832       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1833       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1834       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1835       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1836       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1837      &         'Options in energy minimization:'
1838       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1839      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1840      & 'MinMin:',MinMin,' MinFun:',MinFun,
1841      & ' TolF:',TolF,' RTolF:',RTolF
1842       return
1843       end
1844 c----------------------------------------------------------------------------
1845       subroutine read_angles(kanal,*)
1846       implicit none
1847       include 'DIMENSIONS'
1848       include 'COMMON.GEO'
1849       include 'COMMON.VAR'
1850       include 'COMMON.CHAIN'
1851       include 'COMMON.IOUNITS'
1852       include 'COMMON.CONTROL'
1853       integer i,kanal
1854 c Read angles from input
1855 c
1856        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1857        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1858        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1859        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1860
1861        do i=1,nres
1862 c 9/7/01 avoid 180 deg valence angle
1863         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1864 c
1865         theta(i)=deg2rad*theta(i)
1866         phi(i)=deg2rad*phi(i)
1867         alph(i)=deg2rad*alph(i)
1868         omeg(i)=deg2rad*omeg(i)
1869        enddo
1870       return
1871    10 return1
1872       end
1873 c----------------------------------------------------------------------------
1874       subroutine reada(rekord,lancuch,wartosc,default)
1875       implicit none
1876       character*(*) rekord,lancuch
1877       double precision wartosc,default
1878       integer ilen,iread
1879       external ilen
1880       iread=index(rekord,lancuch)
1881       if (iread.eq.0) then
1882         wartosc=default
1883         return
1884       endif
1885       iread=iread+ilen(lancuch)+1
1886       read (rekord(iread:),*,err=10,end=10) wartosc
1887       return
1888   10  wartosc=default
1889       return
1890       end
1891 c----------------------------------------------------------------------------
1892       subroutine readi(rekord,lancuch,wartosc,default)
1893       implicit none
1894       character*(*) rekord,lancuch
1895       integer wartosc,default
1896       integer ilen,iread
1897       external ilen
1898       iread=index(rekord,lancuch)
1899       if (iread.eq.0) then
1900         wartosc=default
1901         return
1902       endif
1903       iread=iread+ilen(lancuch)+1
1904       read (rekord(iread:),*,err=10,end=10) wartosc
1905       return
1906   10  wartosc=default
1907       return
1908       end
1909 c----------------------------------------------------------------------------
1910       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1911       implicit none
1912       integer dim,i
1913       integer tablica(dim),default
1914       character*(*) rekord,lancuch
1915       character*80 aux
1916       integer ilen,iread
1917       external ilen
1918       do i=1,dim
1919         tablica(i)=default
1920       enddo
1921       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1922       if (iread.eq.0) return
1923       iread=iread+ilen(lancuch)+1
1924       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1925    10 return
1926       end
1927 c----------------------------------------------------------------------------
1928       subroutine multreada(rekord,lancuch,tablica,dim,default)
1929       implicit none
1930       integer dim,i
1931       double precision tablica(dim),default
1932       character*(*) rekord,lancuch
1933       character*80 aux
1934       integer ilen,iread
1935       external ilen
1936       do i=1,dim
1937         tablica(i)=default
1938       enddo
1939       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1940       if (iread.eq.0) return
1941       iread=iread+ilen(lancuch)+1
1942       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1943    10 return
1944       end
1945 c----------------------------------------------------------------------------
1946       subroutine openunits
1947       implicit none
1948       include 'DIMENSIONS'
1949 #ifdef MPI
1950       include 'mpif.h'
1951       integer ierror
1952       character*16 form,nodename
1953       integer nodelen
1954 #endif
1955       include 'COMMON.SETUP'
1956       include 'COMMON.IOUNITS'
1957       include 'COMMON.MD'
1958       include 'COMMON.CONTROL'
1959       integer lenpre,lenpot,ilen,lentmp,npos
1960       external ilen
1961       character*3 out1file_text,ucase
1962       character*3 ll
1963       external ucase
1964       tmpdir=""
1965       out1file_text="YES"
1966 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1967       call getenv_loc("PREFIX",prefix)
1968       pref_orig = prefix
1969       call getenv_loc("POT",pot)
1970       call getenv_loc("DIRTMP",tmpdir)
1971       call getenv_loc("CURDIR",curdir)
1972       call getenv_loc("OUT1FILE",out1file_text)
1973 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1974       out1file_text=ucase(out1file_text)
1975       if (out1file_text(1:1).eq."Y") then
1976         out1file=.true.
1977       else
1978         out1file=fg_rank.gt.0
1979       endif
1980       lenpre=ilen(prefix)
1981       lenpot=ilen(pot)
1982       lentmp=ilen(tmpdir)
1983       if (lentmp.gt.0) then
1984           write (*,'(80(1h!))')
1985           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1986           write (*,'(80(1h!))')
1987           write (*,*)"All output files will be on node /tmp directory."
1988 #ifdef MPI
1989         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1990         if (me.eq.king) then
1991           write (*,*) "The master node is ",nodename
1992         else if (fg_rank.eq.0) then
1993           write (*,*) "I am the CG slave node ",nodename
1994         else
1995           write (*,*) "I am the FG slave node ",nodename
1996         endif
1997 #endif
1998         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1999         lenpre = lentmp+lenpre+1
2000       endif
2001       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2002 C Get the names and open the input files
2003 #if defined(WINIFL) || defined(WINPGI)
2004       open(1,file=pref_orig(:ilen(pref_orig))//
2005      &  '.inp',status='old',readonly,shared)
2006        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2007 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2008 C Get parameter filenames and open the parameter files.
2009       call getenv_loc('BONDPAR',bondname)
2010       open (ibond,file=bondname,status='old',readonly,shared)
2011       call getenv_loc('THETPAR',thetname)
2012       open (ithep,file=thetname,status='old',readonly,shared)
2013       call getenv_loc('ROTPAR',rotname)
2014       open (irotam,file=rotname,status='old',readonly,shared)
2015       call getenv_loc('TORPAR',torname)
2016       open (itorp,file=torname,status='old',readonly,shared)
2017       call getenv_loc('TORDPAR',tordname)
2018       open (itordp,file=tordname,status='old',readonly,shared)
2019       call getenv_loc('FOURIER',fouriername)
2020       open (ifourier,file=fouriername,status='old',readonly,shared)
2021       call getenv_loc('ELEPAR',elename)
2022       open (ielep,file=elename,status='old',readonly,shared)
2023       call getenv_loc('SIDEPAR',sidename)
2024       open (isidep,file=sidename,status='old',readonly,shared)
2025       call getenv_loc('LIPTRANPAR',liptranname)
2026       open (iliptranpar,file=liptranname,status='old',readonly,shared)
2027 #elif (defined CRAY) || (defined AIX)
2028       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2029      &  action='read')
2030 c      print *,"Processor",myrank," opened file 1"
2031       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2032 c      print *,"Processor",myrank," opened file 9"
2033 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2034 C Get parameter filenames and open the parameter files.
2035       call getenv_loc('BONDPAR',bondname)
2036       open (ibond,file=bondname,status='old',action='read')
2037 c      print *,"Processor",myrank," opened file IBOND"
2038       call getenv_loc('THETPAR',thetname)
2039       open (ithep,file=thetname,status='old',action='read')
2040 #ifndef CRYST_THETA
2041       call getenv_loc('THETPARPDB',thetname_pdb)
2042       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2043 #endif
2044 c      print *,"Processor",myrank," opened file ITHEP"
2045       call getenv_loc('ROTPAR',rotname)
2046       open (irotam,file=rotname,status='old',action='read')
2047 #ifndef CRYST_SC
2048       call getenv_loc('ROTPARPDB',rotname_pdb)
2049       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2050 #endif
2051 c      print *,"Processor",myrank," opened file IROTAM"
2052       call getenv_loc('TORPAR',torname)
2053       open (itorp,file=torname,status='old',action='read')
2054 c      print *,"Processor",myrank," opened file ITORP"
2055       call getenv_loc('TORDPAR',tordname)
2056       open (itordp,file=tordname,status='old',action='read')
2057 c      print *,"Processor",myrank," opened file ITORDP"
2058       call getenv_loc('SCCORPAR',sccorname)
2059       open (isccor,file=sccorname,status='old',action='read')
2060 c      print *,"Processor",myrank," opened file ISCCOR"
2061       call getenv_loc('FOURIER',fouriername)
2062       open (ifourier,file=fouriername,status='old',action='read')
2063 c      print *,"Processor",myrank," opened file IFOURIER"
2064       call getenv_loc('ELEPAR',elename)
2065       open (ielep,file=elename,status='old',action='read')
2066 c      print *,"Processor",myrank," opened file IELEP"
2067       call getenv_loc('SIDEPAR',sidename)
2068       open (isidep,file=sidename,status='old',action='read')
2069       call getenv_loc('LIPTRANPAR',liptranname)
2070       open (iliptranpar,file=liptranname,status='old',action='read')
2071 c      print *,"Processor",myrank," opened file ISIDEP"
2072 c      print *,"Processor",myrank," opened parameter files"
2073 #elif (defined G77)
2074       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2075       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2076 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2077 C Get parameter filenames and open the parameter files.
2078       call getenv_loc('BONDPAR',bondname)
2079       open (ibond,file=bondname,status='old')
2080       call getenv_loc('THETPAR',thetname)
2081       open (ithep,file=thetname,status='old')
2082 #ifndef CRYST_THETA
2083       call getenv_loc('THETPARPDB',thetname_pdb)
2084       open (ithep_pdb,file=thetname_pdb,status='old')
2085 #endif
2086       call getenv_loc('ROTPAR',rotname)
2087       open (irotam,file=rotname,status='old')
2088 #ifndef CRYST_SC
2089       call getenv_loc('ROTPARPDB',rotname_pdb)
2090       open (irotam_pdb,file=rotname_pdb,status='old')
2091 #endif
2092       call getenv_loc('TORPAR',torname)
2093       open (itorp,file=torname,status='old')
2094       call getenv_loc('TORDPAR',tordname)
2095       open (itordp,file=tordname,status='old')
2096       call getenv_loc('SCCORPAR',sccorname)
2097       open (isccor,file=sccorname,status='old')
2098       call getenv_loc('FOURIER',fouriername)
2099       open (ifourier,file=fouriername,status='old')
2100       call getenv_loc('ELEPAR',elename)
2101       open (ielep,file=elename,status='old')
2102       call getenv_loc('SIDEPAR',sidename)
2103       open (isidep,file=sidename,status='old')
2104       call getenv_loc('LIPTRANPAR',liptranname)
2105       open (iliptranpar,file=liptranname,status='old')
2106 #else
2107       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2108      &  readonly)
2109        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2110 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2111 C Get parameter filenames and open the parameter files.
2112       call getenv_loc('BONDPAR',bondname)
2113       open (ibond,file=bondname,status='old',readonly)
2114       call getenv_loc('THETPAR',thetname)
2115       open (ithep,file=thetname,status='old',readonly)
2116       call getenv_loc('ROTPAR',rotname)
2117       open (irotam,file=rotname,status='old',readonly)
2118       call getenv_loc('TORPAR',torname)
2119       open (itorp,file=torname,status='old',readonly)
2120       call getenv_loc('TORDPAR',tordname)
2121       open (itordp,file=tordname,status='old',readonly)
2122       call getenv_loc('SCCORPAR',sccorname)
2123       open (isccor,file=sccorname,status='old',readonly)
2124 #ifndef CRYST_THETA
2125       call getenv_loc('THETPARPDB',thetname_pdb)
2126       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2127 #endif
2128       call getenv_loc('FOURIER',fouriername)
2129       open (ifourier,file=fouriername,status='old',readonly)
2130       call getenv_loc('ELEPAR',elename)
2131       open (ielep,file=elename,status='old',readonly)
2132       call getenv_loc('SIDEPAR',sidename)
2133       open (isidep,file=sidename,status='old',readonly)
2134       call getenv_loc('LIPTRANPAR',liptranname)
2135       open (iliptranpar,file=liptranname,status='old',action='read')
2136 #ifndef CRYST_SC
2137       call getenv_loc('ROTPARPDB',rotname_pdb)
2138       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2139 #endif
2140 #endif
2141 #ifdef TUBE
2142       call getenv_loc('TUBEPAR',tubename)
2143 #if defined(WINIFL) || defined(WINPGI)
2144       open (itube,file=tubename,status='old',readonly,shared)
2145 #elif (defined CRAY)  || (defined AIX)
2146       open (itube,file=tubename,status='old',action='read')
2147 #elif (defined G77)
2148       open (itube,file=tubename,status='old')
2149 #else
2150       open (itube,file=tubename,status='old',readonly)
2151 #endif
2152 #endif
2153 #ifndef OLDSCP
2154 C
2155 C 8/9/01 In the newest version SCp interaction constants are read from a file
2156 C Use -DOLDSCP to use hard-coded constants instead.
2157 C
2158       call getenv_loc('SCPPAR',scpname)
2159 #if defined(WINIFL) || defined(WINPGI)
2160       open (iscpp,file=scpname,status='old',readonly,shared)
2161 #elif (defined CRAY)  || (defined AIX)
2162       open (iscpp,file=scpname,status='old',action='read')
2163 #elif (defined G77)
2164       open (iscpp,file=scpname,status='old')
2165 #else
2166       open (iscpp,file=scpname,status='old',readonly)
2167 #endif
2168 #endif
2169       call getenv_loc('PATTERN',patname)
2170 #if defined(WINIFL) || defined(WINPGI)
2171       open (icbase,file=patname,status='old',readonly,shared)
2172 #elif (defined CRAY)  || (defined AIX)
2173       open (icbase,file=patname,status='old',action='read')
2174 #elif (defined G77)
2175       open (icbase,file=patname,status='old')
2176 #else
2177       open (icbase,file=patname,status='old',readonly)
2178 #endif
2179 #ifdef MPI
2180 C Open output file only for CG processes
2181 c      print *,"Processor",myrank," fg_rank",fg_rank
2182       if (fg_rank.eq.0) then
2183
2184       if (nodes.eq.1) then
2185         npos=3
2186       else
2187         npos = dlog10(dfloat(nodes-1))+1
2188       endif
2189       if (npos.lt.3) npos=3
2190       write (liczba,'(i1)') npos
2191       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2192      &  //')'
2193       write (liczba,form) me
2194       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2195      &  liczba(:ilen(liczba))
2196       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2197      &  //'.int'
2198       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2199      &  //'.pdb'
2200       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2201      &  liczba(:ilen(liczba))//'.mol2'
2202       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2203      &  liczba(:ilen(liczba))//'.stat'
2204       if (lentmp.gt.0)
2205      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2206      &      //liczba(:ilen(liczba))//'.stat')
2207       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2208      &  //'.rst'
2209       if(usampl) then
2210           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2211      & liczba(:ilen(liczba))//'.const'
2212       endif
2213
2214       endif
2215 #else
2216       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2217       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2218       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2219       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2220       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2221       if (lentmp.gt.0)
2222      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2223      &    //'.stat')
2224       rest2name=prefix(:ilen(prefix))//'.rst'
2225       if(usampl) then
2226          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2227       endif
2228 #endif
2229 #if defined(AIX) || defined(PGI) || defined(CRAY)
2230       if (me.eq.king .or. .not. out1file) then
2231          open(iout,file=outname,status='unknown')
2232       else
2233          open(iout,file="/dev/null",status="unknown")
2234       endif
2235 c#define DEBUG
2236 #ifdef DEBUG
2237       if (fg_rank.gt.0) then
2238         write (liczba,'(i3.3)') myrank/nfgtasks
2239         write (ll,'(bz,i3.3)') fg_rank
2240         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2241      &   status='unknown')
2242       endif
2243 #endif
2244 c#undef DEBUG
2245       if(me.eq.king) then
2246        open(igeom,file=intname,status='unknown',position='append')
2247        open(ipdb,file=pdbname,status='unknown')
2248        open(imol2,file=mol2name,status='unknown')
2249        open(istat,file=statname,status='unknown',position='append')
2250       else
2251 c1out       open(iout,file=outname,status='unknown')
2252       endif
2253 #else
2254       if (me.eq.king .or. .not.out1file)
2255      &    open(iout,file=outname,status='unknown')
2256 #ifdef DEBUG
2257       if (fg_rank.gt.0) then
2258         write (liczba,'(i3.3)') myrank/nfgtasks
2259         write (ll,'(bz,i3.3)') fg_rank
2260         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2261      &   status='unknown')
2262       endif
2263 #endif
2264       if(me.eq.king) then
2265        open(igeom,file=intname,status='unknown',access='append')
2266        open(ipdb,file=pdbname,status='unknown')
2267        open(imol2,file=mol2name,status='unknown')
2268        open(istat,file=statname,status='unknown',access='append')
2269       else
2270 c1out       open(iout,file=outname,status='unknown')
2271       endif
2272 #endif
2273       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2274       csa_seed=prefix(:lenpre)//'.CSA.seed'
2275       csa_history=prefix(:lenpre)//'.CSA.history'
2276       csa_bank=prefix(:lenpre)//'.CSA.bank'
2277       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2278       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2279       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2280 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2281       csa_int=prefix(:lenpre)//'.int'
2282       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2283       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2284       csa_in=prefix(:lenpre)//'.CSA.in'
2285 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2286 C Write file names
2287       if (me.eq.king)then
2288       write (iout,'(80(1h-))')
2289       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2290       write (iout,'(80(1h-))')
2291       write (iout,*) "Input file                      : ",
2292      &  pref_orig(:ilen(pref_orig))//'.inp'
2293       write (iout,*) "Output file                     : ",
2294      &  outname(:ilen(outname))
2295       write (iout,*)
2296       write (iout,*) "Sidechain potential file        : ",
2297      &  sidename(:ilen(sidename))
2298 #ifndef OLDSCP
2299       write (iout,*) "SCp potential file              : ",
2300      &  scpname(:ilen(scpname))
2301 #endif
2302       write (iout,*) "Electrostatic potential file    : ",
2303      &  elename(:ilen(elename))
2304       write (iout,*) "Cumulant coefficient file       : ",
2305      &  fouriername(:ilen(fouriername))
2306       write (iout,*) "Torsional parameter file        : ",
2307      &  torname(:ilen(torname))
2308       write (iout,*) "Double torsional parameter file : ",
2309      &  tordname(:ilen(tordname))
2310       write (iout,*) "SCCOR parameter file : ",
2311      &  sccorname(:ilen(sccorname))
2312       write (iout,*) "Bond & inertia constant file    : ",
2313      &  bondname(:ilen(bondname))
2314       write (iout,*) "Bending-torsion parameter file  : ",
2315      &  thetname(:ilen(thetname))
2316       write (iout,*) "Rotamer parameter file          : ",
2317      &  rotname(:ilen(rotname))
2318       write (iout,*) "Thetpdb parameter file          : ",
2319      &  thetname_pdb(:ilen(thetname_pdb))
2320       write (iout,*) "Threading database              : ",
2321      &  patname(:ilen(patname))
2322       if (lentmp.ne.0)
2323      &write (iout,*)" DIRTMP                          : ",
2324      &  tmpdir(:lentmp)
2325       write (iout,'(80(1h-))')
2326       endif
2327       return
2328       end
2329 c----------------------------------------------------------------------------
2330       subroutine card_concat(card)
2331       implicit real*8 (a-h,o-z)
2332       include 'DIMENSIONS'
2333       include 'COMMON.IOUNITS'
2334       character*(*) card
2335       character*80 karta,ucase
2336       external ilen
2337       read (inp,'(a)') karta
2338       karta=ucase(karta)
2339       card=' '
2340       do while (karta(80:80).eq.'&')
2341         card=card(:ilen(card)+1)//karta(:79)
2342         read (inp,'(a)') karta
2343         karta=ucase(karta)
2344       enddo
2345       card=card(:ilen(card)+1)//karta
2346       return
2347       end
2348 c------------------------------------------------------------------------------
2349       subroutine readrst
2350       implicit none
2351       include 'DIMENSIONS'
2352       include 'COMMON.CHAIN'
2353       include 'COMMON.IOUNITS'
2354       include 'COMMON.CONTROL'
2355       include 'COMMON.MD'
2356       include 'COMMON.QRESTR'
2357 #ifdef FIVEDIAG
2358        include 'COMMON.LAGRANGE.5diag'
2359 #else
2360        include 'COMMON.LAGRANGE'
2361 #endif
2362       integer i,j
2363       open(irest2,file=rest2name,status='unknown')
2364       read(irest2,*) totT,EK,potE,totE,t_bath
2365       totTafm=totT
2366       do i=1,2*nres
2367          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2368       enddo
2369       do i=1,2*nres
2370          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2371       enddo
2372       if(usampl) then
2373              read (irest2,*) iset
2374       endif
2375       close(irest2)
2376       return
2377       end
2378 c------------------------------------------------------------------------------
2379       subroutine read_fragments
2380       implicit none
2381       include 'DIMENSIONS'
2382 #ifdef MPI
2383       include 'mpif.h'
2384 #endif
2385       include 'COMMON.SETUP'
2386       include 'COMMON.CHAIN'
2387       include 'COMMON.IOUNITS'
2388       include 'COMMON.MD'
2389       include 'COMMON.QRESTR'
2390       include 'COMMON.CONTROL'
2391       integer i
2392       read(inp,*) nset,nfrag,npair,nfrag_back
2393       loc_qlike=(nfrag_back.lt.0)
2394       nfrag_back=iabs(nfrag_back)
2395       if(me.eq.king.or..not.out1file)
2396      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2397      &  " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2398       do iset=1,nset
2399          read(inp,*) mset(iset)
2400        do i=1,nfrag
2401          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2402      &     qinfrag(i,iset)
2403          if(me.eq.king.or..not.out1file)
2404      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2405      &     ifrag(2,i,iset), qinfrag(i,iset)
2406        enddo
2407        do i=1,npair
2408         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2409      &    qinpair(i,iset)
2410         if(me.eq.king.or..not.out1file)
2411      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2412      &    ipair(2,i,iset), qinpair(i,iset)
2413        enddo
2414        if (loc_qlike) then
2415        do i=1,nfrag_back
2416         read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2417      &     wfrag_back(2,i,iset),qin_back(2,i,iset),
2418      &     wfrag_back(3,i,iset),qin_back(3,i,iset),
2419      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2420         if(me.eq.king.or..not.out1file)
2421      &   write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2422      &   wfrag_back(2,i,iset),qin_back(3,i,iset),
2423      &   wfrag_back(3,i,iset),qin_back(3,i,iset),
2424      &   ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2425        enddo
2426        else
2427        do i=1,nfrag_back
2428         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2429      &     wfrag_back(3,i,iset),
2430      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2431         if(me.eq.king.or..not.out1file)
2432      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2433      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2434        enddo
2435        endif
2436       enddo
2437       return
2438       end
2439 C---------------------------------------------------------------------------
2440       subroutine read_afminp
2441       implicit none
2442       include 'DIMENSIONS'
2443 #ifdef MPI
2444       include 'mpif.h'
2445 #endif
2446       include 'COMMON.SETUP'
2447       include 'COMMON.CONTROL'
2448       include 'COMMON.CHAIN'
2449       include 'COMMON.IOUNITS'
2450       include 'COMMON.SBRIDGE'
2451       character*320 afmcard
2452       integer i
2453 c      print *, "wchodze"
2454       call card_concat(afmcard)
2455       call readi(afmcard,"BEG",afmbeg,0)
2456       call readi(afmcard,"END",afmend,0)
2457       call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2458       call reada(afmcard,"VEL",velAFMconst,0.0d0)
2459       print *,'FORCE=' ,forceAFMconst
2460 CCCC NOW PROPERTIES FOR AFM
2461        distafminit=0.0d0
2462        do i=1,3
2463         distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2464        enddo
2465         distafminit=dsqrt(distafminit)
2466 c        print *,'initdist',distafminit
2467       return
2468       end
2469 c-------------------------------------------------------------------------------
2470       subroutine read_saxs_constr
2471       implicit none
2472       include 'DIMENSIONS'
2473 #ifdef MPI
2474       include 'mpif.h'
2475 #endif
2476       include 'COMMON.SETUP'
2477       include 'COMMON.CONTROL'
2478       include 'COMMON.SAXS'
2479       include 'COMMON.CHAIN'
2480       include 'COMMON.IOUNITS'
2481       include 'COMMON.SBRIDGE'
2482       double precision cm(3),cnorm
2483       integer i,j
2484 c      read(inp,*) nsaxs
2485       write (iout,*) "Calling read_saxs nsaxs",nsaxs
2486       call flush(iout)
2487       if (saxs_mode.eq.0) then
2488 c SAXS distance distribution
2489       do i=1,nsaxs
2490         read(inp,*) distsaxs(i),Psaxs(i)
2491       enddo
2492       Cnorm = 0.0d0
2493       do i=1,nsaxs
2494         Cnorm = Cnorm + Psaxs(i)
2495       enddo
2496       write (iout,*) "Cnorm",Cnorm
2497       do i=1,nsaxs
2498         Psaxs(i)=Psaxs(i)/Cnorm
2499       enddo
2500       write (iout,*) "Normalized distance distribution from SAXS"
2501       do i=1,nsaxs
2502         write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2503       enddo
2504       Wsaxs0=0.0d0
2505       do i=1,nsaxs
2506         Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2507       enddo
2508       write (iout,*) "Wsaxs0",Wsaxs0
2509       else
2510 c SAXS "spheres".
2511       do i=1,nsaxs
2512         read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2513       enddo
2514       do j=1,3
2515         cm(j)=0.0d0
2516       enddo
2517       do i=1,nsaxs
2518         do j=1,3
2519           cm(j)=cm(j)+Csaxs(j,i)
2520         enddo
2521       enddo
2522       do j=1,3
2523         cm(j)=cm(j)/nsaxs
2524       enddo
2525       do i=1,nsaxs
2526         do j=1,3
2527           Csaxs(j,i)=Csaxs(j,i)-cm(j)
2528         enddo
2529       enddo
2530       write (iout,*) "SAXS sphere coordinates"
2531       do i=1,nsaxs
2532         write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2533       enddo
2534       endif
2535       return
2536       end
2537
2538 c-------------------------------------------------------------------------------
2539       subroutine read_dist_constr
2540       implicit none
2541       include 'DIMENSIONS'
2542 #ifdef MPI
2543       include 'mpif.h'
2544 #endif
2545       include 'COMMON.SETUP'
2546       include 'COMMON.CONTROL'
2547       include 'COMMON.CHAIN'
2548       include 'COMMON.IOUNITS'
2549       include 'COMMON.SBRIDGE'
2550       include 'COMMON.INTERACT'
2551       integer i,j,k,ii,jj,itemp,link_type,iiend,jjend,kk
2552       integer nfrag_,npair_,ndist_,ifrag_(2,100),ipair_(2,1000)
2553       double precision wfrag_(100),wpair_(1000)
2554       double precision ddjk,dist,dist_cut,fordepthmax
2555       character*5000 controlcard
2556       logical normalize,next
2557       integer restr_type
2558       double precision scal_bfac
2559       double precision xlink(4,0:4) /
2560 c           a          b       c     sigma
2561      &   0.0d0,0.0d0,0.0d0,0.0d0,                             ! default, no xlink potential
2562      &   0.00305218d0,9.46638d0,4.68901d0,4.74347d0,          ! ZL
2563      &   0.00214928d0,12.7517d0,0.00375009d0,6.13477d0,       ! ADH
2564      &   0.00184547d0,11.2678d0,0.00140292d0,7.00868d0,       ! PDH
2565      &   0.000161786d0,6.29273d0,4.40993d0,7.13956d0    /     ! DSS
2566 c      print *, "WCHODZE"
2567       write (iout,*) "Calling read_dist_constr"
2568 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2569 c      call flush(iout)
2570       restr_on_coord=.false.
2571       next=.true.
2572
2573       npeak=0
2574       ipeak=0
2575       nhpb_peak=0
2576
2577       DO WHILE (next)
2578
2579       call card_concat(controlcard)
2580       next = index(controlcard,"NEXT").gt.0
2581       call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2582       write (iout,*) "restr_type",restr_type
2583       call readi(controlcard,"NFRAG",nfrag_,0)
2584       call readi(controlcard,"NFRAG",nfrag_,0)
2585       call readi(controlcard,"NPAIR",npair_,0)
2586       call readi(controlcard,"NDIST",ndist_,0)
2587       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2588       call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2589       if (restr_type.eq.10)
2590      &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2591       if (restr_type.eq.12)
2592      &  call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2593       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2594       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2595       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2596       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2597       normalize = index(controlcard,"NORMALIZE").gt.0
2598       write (iout,*) "WBOLTZD",wboltzd
2599       write (iout,*) "SCAL_PEAK",scal_peak
2600       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2601       write (iout,*) "IFRAG"
2602       do i=1,nfrag_
2603         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2604       enddo
2605       write (iout,*) "IPAIR"
2606       do i=1,npair_
2607         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2608       enddo
2609       if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2610      & write (iout,*)
2611      &   "Distance restraints as generated from reference structure"
2612       do i=1,nfrag_
2613         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2614         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2615      &    ifrag_(2,i)=nstart_sup+nsup-1
2616 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2617 c        call flush(iout)
2618         if (wfrag_(i).eq.0.0d0) cycle
2619         do j=ifrag_(1,i),ifrag_(2,i)-1
2620           do k=j+1,ifrag_(2,i)
2621 c            write (iout,*) "j",j," k",k
2622             ddjk=dist(j,k)
2623             if (restr_type.eq.1) then
2624               nhpb=nhpb+1
2625               irestr_type(nhpb)=1
2626               ihpb(nhpb)=j
2627               jhpb(nhpb)=k
2628               dhpb(nhpb)=ddjk
2629               forcon(nhpb)=wfrag_(i)
2630             else if (constr_dist.eq.2) then
2631               if (ddjk.le.dist_cut) then
2632                 nhpb=nhpb+1
2633                 irestr_type(nhpb)=1
2634                 ihpb(nhpb)=j
2635                 jhpb(nhpb)=k
2636                 dhpb(nhpb)=ddjk
2637                 forcon(nhpb)=wfrag_(i)
2638               endif
2639             else if (restr_type.eq.3) then
2640               nhpb=nhpb+1
2641               irestr_type(nhpb)=1
2642               ihpb(nhpb)=j
2643               jhpb(nhpb)=k
2644               dhpb(nhpb)=ddjk
2645               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2646             endif
2647 #ifdef MPI
2648             if (.not.out1file .or. me.eq.king)
2649      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2650      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2651 #else
2652             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2653      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2654 #endif
2655           enddo
2656         enddo
2657       enddo
2658       do i=1,npair_
2659         if (wpair_(i).eq.0.0d0) cycle
2660         ii = ipair_(1,i)
2661         jj = ipair_(2,i)
2662         if (ii.gt.jj) then
2663           itemp=ii
2664           ii=jj
2665           jj=itemp
2666         endif
2667         do j=ifrag_(1,ii),ifrag_(2,ii)
2668           do k=ifrag_(1,jj),ifrag_(2,jj)
2669             ddjk=dist(j,k)
2670             if (restr_type.eq.1) then
2671               nhpb=nhpb+1
2672               irestr_type(nhpb)=1
2673               ihpb(nhpb)=j
2674               jhpb(nhpb)=k
2675               dhpb(nhpb)=ddjk
2676               forcon(nhpb)=wpair_(i)
2677             else if (constr_dist.eq.2) then
2678               if (ddjk.le.dist_cut) then
2679                 nhpb=nhpb+1
2680                 irestr_type(nhpb)=1
2681                 ihpb(nhpb)=j
2682                 jhpb(nhpb)=k
2683                 dhpb(nhpb)=ddjk
2684                 forcon(nhpb)=wpair_(i)
2685               endif
2686             else if (restr_type.eq.3) then
2687               nhpb=nhpb+1
2688               irestr_type(nhpb)=1
2689               ihpb(nhpb)=j
2690               jhpb(nhpb)=k
2691               dhpb(nhpb)=ddjk
2692               forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2693             endif
2694 #ifdef MPI
2695             if (.not.out1file .or. me.eq.king)
2696      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2697      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2698 #else
2699             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2700      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2701 #endif
2702           enddo
2703         enddo
2704       enddo
2705
2706 c      print *,ndist_
2707       write (iout,*) "Distance restraints as read from input"
2708       do i=1,ndist_
2709         if (restr_type.eq.12) then
2710           read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2711      &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2712      &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2713      &    fordepth_peak(nhpb_peak+1),npeak
2714 c          write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2715 c     &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2716 c     &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2717 c     &    fordepth_peak(nhpb_peak+1),npeak
2718           if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2719      &      fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2720           nhpb_peak=nhpb_peak+1
2721           irestr_type_peak(nhpb_peak)=12
2722           if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2723           ipeak(2,npeak)=i
2724 #ifdef MPI
2725           if (.not.out1file .or. me.eq.king)
2726      &    write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2727      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2728      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2729      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2730      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2731 #else
2732           write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2733      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2734      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2735      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2736      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2737 #endif
2738           if (ibecarb_peak(nhpb_peak).eq.3) then
2739             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2740           else if (ibecarb_peak(nhpb_peak).eq.2) then
2741             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2742           else if (ibecarb_peak(nhpb_peak).eq.1) then
2743             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2744             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2745           endif
2746         else if (restr_type.eq.11) then
2747           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2748      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2749 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2750           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2751           nhpb=nhpb+1
2752           irestr_type(nhpb)=11
2753 #ifdef MPI
2754           if (.not.out1file .or. me.eq.king)
2755      &    write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2756      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2757      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2758 #else
2759           write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2760      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2761      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2762 #endif
2763 c          if (ibecarb(nhpb).gt.0) then
2764 c            ihpb(nhpb)=ihpb(nhpb)+nres
2765 c            jhpb(nhpb)=jhpb(nhpb)+nres
2766 c          endif
2767          if (ibecarb(nhpb).eq.3) then
2768             ihpb(nhpb)=ihpb(nhpb)+nres
2769           else if (ibecarb(nhpb).eq.2) then
2770             ihpb(nhpb)=ihpb(nhpb)+nres
2771           else if (ibecarb(nhpb).eq.1) then
2772             ihpb(nhpb)=ihpb(nhpb)+nres
2773             jhpb(nhpb)=jhpb(nhpb)+nres
2774           endif
2775         else if (restr_type.eq.10) then
2776 c Cross-lonk Markov-like potential
2777           call card_concat(controlcard)
2778           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2779           call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2780           ibecarb(nhpb+1)=0
2781           if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2782           if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2783           if (index(controlcard,"ZL").gt.0) then
2784             link_type=1
2785           else if (index(controlcard,"ADH").gt.0) then
2786             link_type=2
2787           else if (index(controlcard,"PDH").gt.0) then
2788             link_type=3
2789           else if (index(controlcard,"DSS").gt.0) then
2790             link_type=4
2791           else
2792             link_type=0
2793           endif
2794           call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2795      &       xlink(1,link_type))
2796           call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2797      &       xlink(2,link_type))
2798           call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2799      &       xlink(3,link_type))
2800           call reada(controlcard,"SIGMA",forcon(nhpb+1),
2801      &       xlink(4,link_type))
2802           call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2803 c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2804 c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2805           if (forcon(nhpb+1).le.0.0d0 .or.
2806      &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2807           nhpb=nhpb+1
2808           irestr_type(nhpb)=10
2809           if (ibecarb(nhpb).eq.3) then
2810             jhpb(nhpb)=jhpb(nhpb)+nres
2811           else if (ibecarb(nhpb).eq.2) then
2812             ihpb(nhpb)=ihpb(nhpb)+nres
2813           else if (ibecarb(nhpb).eq.1) then
2814             ihpb(nhpb)=ihpb(nhpb)+nres
2815             jhpb(nhpb)=jhpb(nhpb)+nres
2816           endif
2817 #ifdef MPI
2818           if (.not.out1file .or. me.eq.king)
2819      &    write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2820      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2821      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2822      &     irestr_type(nhpb)
2823 #else
2824           write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2825      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2826      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2827      &     irestr_type(nhpb)
2828 #endif
2829         else
2830 C        print *,"in else"
2831           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2832      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2833           if (forcon(nhpb+1).gt.0.0d0) then
2834           nhpb=nhpb+1
2835           if (dhpb1(nhpb).eq.0.0d0) then
2836             irestr_type(nhpb)=1
2837           else
2838             irestr_type(nhpb)=2
2839           endif
2840           if (ibecarb(nhpb).eq.3) then
2841             jhpb(nhpb)=jhpb(nhpb)+nres
2842           else if (ibecarb(nhpb).eq.2) then
2843             ihpb(nhpb)=ihpb(nhpb)+nres
2844           else if (ibecarb(nhpb).eq.1) then
2845             ihpb(nhpb)=ihpb(nhpb)+nres
2846             jhpb(nhpb)=jhpb(nhpb)+nres
2847           endif
2848           if (dhpb(nhpb).eq.0.0d0)
2849      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2850         endif
2851 #ifdef MPI
2852           if (.not.out1file .or. me.eq.king)
2853      &    write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2854      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2855 #else
2856           write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2857      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2858 #endif
2859         endif
2860 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2861 C        if (forcon(nhpb+1).gt.0.0d0) then
2862 C          nhpb=nhpb+1
2863 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2864       enddo
2865
2866       if (restr_type.eq.4) then
2867         write (iout,*) "The BFAC array"
2868         do i=nnt,nct
2869           write (iout,'(i5,f10.5)') i,bfac(i)
2870         enddo
2871         do i=nnt,nct
2872           if (itype(i).eq.ntyp1) cycle
2873           do j=nnt,i-1
2874             if (itype(j).eq.ntyp1) cycle
2875             if (itype(i).eq.10) then
2876               iiend=0
2877             else
2878               iiend=1
2879             endif
2880             if (itype(j).eq.10) then
2881               jjend=0
2882             else
2883               jjend=1
2884             endif
2885             kk=0
2886             do ii=0,iiend
2887             do jj=0,jjend
2888             nhpb=nhpb+1
2889             irestr_type(nhpb)=1
2890             forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2891             irestr_type(nhpb)=1
2892             ibecarb(nhpb)=kk
2893             if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2894             ihpb(nhpb)=i+nres*ii
2895             jhpb(nhpb)=j+nres*jj
2896             dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2897 #ifdef MPI
2898             if (.not.out1file .or. me.eq.king) then
2899             write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2900      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2901      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2902      &       irestr_type(nhpb)
2903             endif
2904 #else
2905             write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2906      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2907      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2908      &       irestr_type(nhpb)
2909 #endif
2910             kk=kk+1
2911           enddo
2912           enddo
2913           enddo
2914         enddo
2915       endif
2916
2917       if (restr_type.eq.5) then
2918         restr_on_coord=.true.
2919         do i=nnt,nct
2920           if (itype(i).eq.ntyp1) cycle
2921           bfac(i)=(scal_bfac/bfac(i))**2
2922         enddo
2923       endif
2924
2925       ENDDO ! next
2926
2927       fordepthmax=0.0d0
2928       if (normalize) then
2929         do i=nss+1,nhpb
2930           if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
2931      &      fordepthmax=fordepth(i)
2932         enddo
2933         do i=nss+1,nhpb
2934           if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
2935         enddo
2936       endif
2937       return
2938       end
2939 c-------------------------------------------------------------------------------
2940       subroutine read_constr_homology
2941       implicit none
2942       include 'DIMENSIONS'
2943 #ifdef MPI
2944       include 'mpif.h'
2945 #endif
2946       include 'COMMON.SETUP'
2947       include 'COMMON.CONTROL'
2948       include 'COMMON.HOMOLOGY'
2949       include 'COMMON.CHAIN'
2950       include 'COMMON.IOUNITS'
2951       include 'COMMON.MD'
2952       include 'COMMON.QRESTR'
2953       include 'COMMON.GEO'
2954       include 'COMMON.INTERACT'
2955       include 'COMMON.NAMES'
2956 c
2957 c For new homol impl
2958 c
2959       include 'COMMON.VAR'
2960 c
2961
2962 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2963 c    &                 dist_cut
2964 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2965 c    &    sigma_odl_temp(maxres,maxres,max_template)
2966       character*2 kic2
2967       character*24 model_ki_dist, model_ki_angle
2968       character*500 controlcard
2969       integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
2970      & ik,iistart,iishift
2971       integer ilen
2972       external ilen
2973       logical liiflag
2974 c
2975 c     FP - Nov. 2014 Temporary specifications for new vars
2976 c
2977       double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
2978      &                 rescore3_tmp
2979       double precision, dimension (max_template,maxres) :: rescore
2980       double precision, dimension (max_template,maxres) :: rescore2
2981       double precision, dimension (max_template,maxres) :: rescore3
2982       double precision distal
2983       character*24 pdbfile,tpl_k_rescore
2984 c -----------------------------------------------------------------
2985 c Reading multiple PDB ref structures and calculation of retraints
2986 c not using pre-computed ones stored in files model_ki_{dist,angle}
2987 c FP (Nov., 2014)
2988 c -----------------------------------------------------------------
2989 c
2990 c
2991 c Alternative: reading from input
2992       call card_concat(controlcard)
2993       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2994       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2995       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2996       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2997       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2998       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
2999       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3000       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3001       start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3002       if(.not.read2sigma.and.start_from_model) then
3003           if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3004      &      write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3005           start_from_model=.false.
3006       endif
3007       if(start_from_model .and. (me.eq.king .or. .not. out1file))
3008      &    write(iout,*) 'START_FROM_MODELS is ON'
3009       if(start_from_model .and. rest) then
3010         if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3011          write(iout,*) 'START_FROM_MODELS is OFF'
3012          write(iout,*) 'remove restart keyword from input'
3013         endif
3014       endif
3015       if (homol_nset.gt.1)then
3016          call card_concat(controlcard)
3017          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3018          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3019           write(iout,*) "iset homology_weight "
3020           do i=1,homol_nset
3021            write(iout,*) i,waga_homology(i)
3022           enddo
3023          endif
3024          iset=mod(kolor,homol_nset)+1
3025       else
3026        iset=1
3027        waga_homology(1)=1.0
3028       endif
3029
3030 cd      write (iout,*) "nnt",nnt," nct",nct
3031 cd      call flush(iout)
3032
3033
3034       lim_odl=0
3035       lim_dih=0
3036 c
3037 c      write(iout,*) 'nnt=',nnt,'nct=',nct
3038 c
3039       do i = nnt,nct
3040         do k=1,constr_homology
3041          idomain(k,i)=0
3042         enddo
3043       enddo
3044
3045       ii=0
3046       do i = nnt,nct-2
3047         do j=i+2,nct
3048         ii=ii+1
3049         ii_in_use(ii)=0
3050         enddo
3051       enddo
3052
3053       if (read_homol_frag) then
3054         call read_klapaucjusz
3055       else
3056
3057       do k=1,constr_homology
3058
3059         read(inp,'(a)') pdbfile
3060         if(me.eq.king .or. .not. out1file)
3061      &    write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3062      &  pdbfile(:ilen(pdbfile))
3063         open(ipdbin,file=pdbfile,status='old',err=33)
3064         goto 34
3065   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
3066      &  pdbfile(:ilen(pdbfile))
3067         stop
3068   34    continue
3069 c        print *,'Begin reading pdb data'
3070 c
3071 c Files containing res sim or local scores (former containing sigmas)
3072 c
3073
3074         write(kic2,'(bz,i2.2)') k
3075
3076         tpl_k_rescore="template"//kic2//".sco"
3077
3078         unres_pdb=.false.
3079         if (read2sigma) then
3080           call readpdb_template(k)
3081         else
3082           call readpdb
3083         endif
3084 c
3085 c     Distance restraints
3086 c
3087 c          ... --> odl(k,ii)
3088 C Copy the coordinates from reference coordinates (?)
3089         do i=1,2*nres
3090           do j=1,3
3091             c(j,i)=cref(j,i)
3092 c           write (iout,*) "c(",j,i,") =",c(j,i)
3093           enddo
3094         enddo
3095 c
3096 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3097 c
3098 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3099           open (ientin,file=tpl_k_rescore,status='old')
3100           if (nnt.gt.1) rescore(k,1)=0.0d0
3101           do irec=nnt,nct ! loop for reading res sim
3102             if (read2sigma) then
3103              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3104      &                                rescore3_tmp,idomain_tmp
3105              i_tmp=i_tmp+nnt-1
3106              idomain(k,i_tmp)=idomain_tmp
3107              rescore(k,i_tmp)=rescore_tmp
3108              rescore2(k,i_tmp)=rescore2_tmp
3109              rescore3(k,i_tmp)=rescore3_tmp
3110              if (.not. out1file .or. me.eq.king)
3111      &        write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3112      &                      i_tmp,rescore2_tmp,rescore_tmp,
3113      &                                rescore3_tmp,idomain_tmp
3114             else
3115              idomain(k,irec)=1
3116              read (ientin,*,end=1401) rescore_tmp
3117
3118 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3119              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3120 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3121             endif
3122           enddo
3123  1401   continue
3124         close (ientin)
3125         if (waga_dist.ne.0.0d0) then
3126           ii=0
3127           do i = nnt,nct-2
3128             do j=i+2,nct
3129
3130               x12=c(1,i)-c(1,j)
3131               y12=c(2,i)-c(2,j)
3132               z12=c(3,i)-c(3,j)
3133               distal=dsqrt(x12*x12+y12*y12+z12*z12)
3134 c              write (iout,*) k,i,j,distal,dist2_cut
3135
3136             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3137      &            .and. distal.le.dist2_cut ) then
3138
3139               ii=ii+1
3140               ii_in_use(ii)=1
3141               l_homo(k,ii)=.true.
3142
3143 c             write (iout,*) "k",k
3144 c             write (iout,*) "i",i," j",j," constr_homology",
3145 c    &                       constr_homology
3146               ires_homo(ii)=i
3147               jres_homo(ii)=j
3148               odl(k,ii)=distal
3149               if (read2sigma) then
3150                 sigma_odl(k,ii)=0
3151                 do ik=i,j
3152                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3153                 enddo
3154                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3155                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3156      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3157               else
3158                 if (odl(k,ii).le.dist_cut) then
3159                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3160                 else
3161 #ifdef OLDSIGMA
3162                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3163      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3164 #else
3165                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3166      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3167 #endif
3168                 endif
3169               endif
3170               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3171             else
3172               ii=ii+1
3173               l_homo(k,ii)=.false.
3174             endif
3175             enddo
3176           enddo
3177         lim_odl=ii
3178         endif
3179 c
3180 c     Theta, dihedral and SC retraints
3181 c
3182         if (waga_angle.gt.0.0d0) then
3183 c         open (ientin,file=tpl_k_sigma_dih,status='old')
3184 c         do irec=1,maxres-3 ! loop for reading sigma_dih
3185 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3186 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3187 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3188 c    &                            sigma_dih(k,i+nnt-1)
3189 c         enddo
3190 c1402   continue
3191 c         close (ientin)
3192           do i = nnt+3,nct
3193             if (idomain(k,i).eq.0) then
3194                sigma_dih(k,i)=0.0
3195                cycle
3196             endif
3197             dih(k,i)=phiref(i) ! right?
3198 c           read (ientin,*) sigma_dih(k,i) ! original variant
3199 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
3200 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3201 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
3202 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
3203 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
3204
3205             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3206      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
3207 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3208 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
3209 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3210 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3211 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
3212 c   Instead of res sim other local measure of b/b str reliability possible
3213             if (sigma_dih(k,i).ne.0)
3214      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3215 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3216           enddo
3217           lim_dih=nct-nnt-2
3218         endif
3219
3220         if (waga_theta.gt.0.0d0) then
3221 c         open (ientin,file=tpl_k_sigma_theta,status='old')
3222 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3223 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3224 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3225 c    &                              sigma_theta(k,i+nnt-1)
3226 c         enddo
3227 c1403   continue
3228 c         close (ientin)
3229
3230           do i = nnt+2,nct ! right? without parallel.
3231 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
3232 c         do i=ithet_start,ithet_end ! with FG parallel.
3233              if (idomain(k,i).eq.0) then
3234               sigma_theta(k,i)=0.0
3235               cycle
3236              endif
3237              thetatpl(k,i)=thetaref(i)
3238 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3239 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
3240 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
3241 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
3242 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
3243              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3244      &                        rescore(k,i-2))/3.0
3245 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3246              if (sigma_theta(k,i).ne.0)
3247      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3248
3249 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3250 c                             rescore(k,i-2) !  right expression ?
3251 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3252           enddo
3253         endif
3254
3255         if (waga_d.gt.0.0d0) then
3256 c       open (ientin,file=tpl_k_sigma_d,status='old')
3257 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3258 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3259 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3260 c    &                          sigma_d(k,i+nnt-1)
3261 c         enddo
3262 c1404   continue
3263
3264           do i = nnt,nct ! right? without parallel.
3265 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
3266 c         do i=loc_start,loc_end ! with FG parallel.
3267                if (itype(i).eq.10) cycle
3268                if (idomain(k,i).eq.0 ) then
3269                   sigma_d(k,i)=0.0
3270                   cycle
3271                endif
3272                xxtpl(k,i)=xxref(i)
3273                yytpl(k,i)=yyref(i)
3274                zztpl(k,i)=zzref(i)
3275 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3276 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3277 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3278 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
3279                sigma_d(k,i)=rescore3(k,i) !  right expression ?
3280                if (sigma_d(k,i).ne.0)
3281      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3282
3283 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
3284 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3285 c              read (ientin,*) sigma_d(k,i) ! 1st variant
3286           enddo
3287         endif
3288       enddo
3289 c
3290 c remove distance restraints not used in any model from the list
3291 c shift data in all arrays
3292 c
3293       if (waga_dist.ne.0.0d0) then
3294         ii=0
3295         liiflag=.true.
3296         do i=nnt,nct-2
3297          do j=i+2,nct
3298           ii=ii+1
3299           if (ii_in_use(ii).eq.0.and.liiflag) then
3300             liiflag=.false.
3301             iistart=ii
3302           endif
3303           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3304      &                   .not.liiflag.and.ii.eq.lim_odl) then
3305              if (ii.eq.lim_odl) then
3306               iishift=ii-iistart+1
3307              else
3308               iishift=ii-iistart
3309              endif
3310              liiflag=.true.
3311              do ki=iistart,lim_odl-iishift
3312               ires_homo(ki)=ires_homo(ki+iishift)
3313               jres_homo(ki)=jres_homo(ki+iishift)
3314               ii_in_use(ki)=ii_in_use(ki+iishift)
3315               do k=1,constr_homology
3316                odl(k,ki)=odl(k,ki+iishift)
3317                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3318                l_homo(k,ki)=l_homo(k,ki+iishift)
3319               enddo
3320              enddo
3321              ii=ii-iishift
3322              lim_odl=lim_odl-iishift
3323           endif
3324          enddo
3325         enddo
3326       endif
3327
3328       endif ! .not. klapaucjusz
3329
3330       if (constr_homology.gt.0) call homology_partition
3331       if (constr_homology.gt.0) call init_int_table
3332 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3333 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3334 c
3335 c Print restraints
3336 c
3337       if (.not.out_template_restr) return
3338 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3339       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3340        write (iout,*) "Distance restraints from templates"
3341        do ii=1,lim_odl
3342        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3343      &  ii,ires_homo(ii),jres_homo(ii),
3344      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3345      &  ki=1,constr_homology)
3346        enddo
3347        write (iout,*) "Dihedral angle restraints from templates"
3348        do i=nnt+3,nct
3349         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3350      &      (rad2deg*dih(ki,i),
3351      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3352        enddo
3353        write (iout,*) "Virtual-bond angle restraints from templates"
3354        do i=nnt+2,nct
3355         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3356      &      (rad2deg*thetatpl(ki,i),
3357      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3358        enddo
3359        write (iout,*) "SC restraints from templates"
3360        do i=nnt,nct
3361         write(iout,'(i5,100(4f8.2,4x))') i,
3362      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3363      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3364        enddo
3365       endif
3366 c -----------------------------------------------------------------
3367       return
3368       end
3369 c----------------------------------------------------------------------
3370 #ifdef WINIFL
3371       subroutine flush(iu)
3372       return
3373       end
3374 #endif
3375 #ifdef AIX
3376       subroutine flush(iu)
3377       call flush_(iu)
3378       return
3379       end
3380 #endif
3381 c------------------------------------------------------------------------------
3382       subroutine copy_to_tmp(source)
3383       implicit none
3384       include "DIMENSIONS"
3385       include "COMMON.IOUNITS"
3386       character*(*) source
3387       character* 256 tmpfile
3388       integer ilen
3389       external ilen
3390       logical ex
3391       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3392       inquire(file=tmpfile,exist=ex)
3393       if (ex) then
3394         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3395      &   " to temporary directory..."
3396         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3397         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3398       endif
3399       return
3400       end
3401 c------------------------------------------------------------------------------
3402       subroutine move_from_tmp(source)
3403       implicit none
3404       include "DIMENSIONS"
3405       include "COMMON.IOUNITS"
3406       character*(*) source
3407       integer ilen
3408       external ilen
3409       write (*,*) "Moving ",source(:ilen(source)),
3410      & " from temporary directory to working directory"
3411       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3412       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3413       return
3414       end
3415 c------------------------------------------------------------------------------
3416       subroutine random_init(seed)
3417 C
3418 C Initialize random number generator
3419 C
3420       implicit none
3421       include 'DIMENSIONS'
3422 #ifdef MPI
3423       include 'mpif.h'
3424       logical OKRandom, prng_restart
3425       real*8  r1
3426       integer iseed_array(4)
3427       integer error_msg,ierr
3428 #endif
3429       include 'COMMON.IOUNITS'
3430       include 'COMMON.TIME1'
3431       include 'COMMON.THREAD'
3432       include 'COMMON.SBRIDGE'
3433       include 'COMMON.CONTROL'
3434       include 'COMMON.MCM'
3435       include 'COMMON.MAP'
3436       include 'COMMON.HEADER'
3437       include 'COMMON.CSA'
3438       include 'COMMON.CHAIN'
3439       include 'COMMON.MUCA'
3440       include 'COMMON.MD'
3441       include 'COMMON.FFIELD'
3442       include 'COMMON.SETUP'
3443       integer i,iseed
3444       double precision seed,ran_number
3445       iseed=-dint(dabs(seed))
3446       if (iseed.eq.0) then
3447         write (iout,'(/80(1h*)/20x,a/80(1h*))')
3448      &    'Random seed undefined. The program will stop.'
3449         write (*,'(/80(1h*)/20x,a/80(1h*))')
3450      &    'Random seed undefined. The program will stop.'
3451 #ifdef MPI
3452         call mpi_finalize(mpi_comm_world,ierr)
3453 #endif
3454         stop 'Bad random seed.'
3455       endif
3456 #ifdef MPI
3457       if (fg_rank.eq.0) then
3458       seed=seed*(me+1)+1
3459 #ifdef AMD64
3460       if(me.eq.king)
3461      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3462       OKRandom = prng_restart(me,iseed)
3463 #else
3464       do i=1,4
3465        tmp=65536.0d0**(4-i)
3466        iseed_array(i) = dint(seed/tmp)
3467        seed=seed-iseed_array(i)*tmp
3468       enddo
3469       if(me.eq.king)
3470      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3471      &                 (iseed_array(i),i=1,4)
3472       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3473      &                 (iseed_array(i),i=1,4)
3474       OKRandom = prng_restart(me,iseed_array)
3475 #endif
3476       if (OKRandom) then
3477 c        r1 = prng_next(me)
3478         r1=ran_number(0.0D0,1.0D0)
3479         if(me.eq.king)
3480      &   write (iout,*) 'ran_num',r1
3481         if (r1.lt.0.0d0) OKRandom=.false.
3482       endif
3483       if (.not.OKRandom) then
3484         write (iout,*) 'PRNG IS NOT WORKING!!!'
3485         print *,'PRNG IS NOT WORKING!!!'
3486         if (me.eq.0) then
3487          call flush(iout)
3488          call mpi_abort(mpi_comm_world,error_msg,ierr)
3489          stop
3490         else
3491          write (iout,*) 'too many processors for parallel prng'
3492          write (*,*) 'too many processors for parallel prng'
3493          call flush(iout)
3494          stop
3495         endif
3496       endif
3497       endif
3498 #else
3499       call vrndst(iseed)
3500       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3501 #endif
3502       return
3503       end
3504 c----------------------------------------------------------------------
3505       subroutine read_klapaucjusz
3506       implicit none
3507       include 'DIMENSIONS'
3508 #ifdef MPI
3509       include 'mpif.h'
3510 #endif
3511       include 'COMMON.SETUP'
3512       include 'COMMON.CONTROL'
3513       include 'COMMON.HOMOLOGY'
3514       include 'COMMON.CHAIN'
3515       include 'COMMON.IOUNITS'
3516       include 'COMMON.MD'
3517       include 'COMMON.GEO'
3518       include 'COMMON.INTERACT'
3519       include 'COMMON.NAMES'
3520       character*256 fragfile
3521       integer ninclust(maxclust),inclust(max_template,maxclust),
3522      &  nresclust(maxclust),iresclust(maxres,maxclust),nclust
3523
3524       character*2 kic2
3525       character*24 model_ki_dist, model_ki_angle
3526       character*500 controlcard
3527       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp,
3528      & ik,ll,ii,kk,iistart,iishift,lim_xx
3529       double precision distal
3530       logical lprn /.true./
3531       integer ilen
3532       external ilen
3533       logical liiflag
3534 c
3535 c
3536       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3537       double precision, dimension (max_template,maxres) :: rescore
3538       double precision, dimension (max_template,maxres) :: rescore2
3539       character*24 pdbfile,tpl_k_rescore
3540
3541 c
3542 c For new homol impl
3543 c
3544       include 'COMMON.VAR'
3545 c
3546       call getenv("FRAGFILE",fragfile)
3547       open(ientin,file=fragfile,status="old",err=10)
3548       read(ientin,*) constr_homology,nclust
3549       l_homo = .false.
3550       sigma_theta=0.0
3551       sigma_d=0.0
3552       sigma_dih=0.0
3553 c Read pdb files
3554       do k=1,constr_homology
3555         read(ientin,'(a)') pdbfile
3556         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3557      &  pdbfile(:ilen(pdbfile))
3558         open(ipdbin,file=pdbfile,status='old',err=33)
3559         goto 34
3560   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
3561      &  pdbfile(:ilen(pdbfile))
3562         stop
3563   34    continue
3564         unres_pdb=.false.
3565         call readpdb_template(k)
3566         do i=1,nres
3567           rescore(k,i)=0.2d0
3568           rescore2(k,i)=1.0d0
3569         enddo
3570       enddo
3571 c Read clusters
3572       do i=1,nclust
3573         read(ientin,*) ninclust(i),nresclust(i)
3574         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3575         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3576       enddo
3577 c
3578 c Loop over clusters
3579 c
3580       do l=1,nclust
3581         do ll = 1,ninclust(l)
3582
3583         k = inclust(ll,l)
3584         do i=1,nres
3585           idomain(k,i)=0
3586         enddo
3587         do i=1,nresclust(l)
3588           if (nnt.gt.1)  then
3589             idomain(k,iresclust(i,l)+1) = 1
3590           else
3591             idomain(k,iresclust(i,l)) = 1
3592           endif
3593         enddo
3594 c
3595 c     Distance restraints
3596 c
3597 c          ... --> odl(k,ii)
3598 C Copy the coordinates from reference coordinates (?)
3599         do i=1,2*nres
3600           do j=1,3
3601             c(j,i)=chomo(j,i,k)
3602 c           write (iout,*) "c(",j,i,") =",c(j,i)
3603           enddo
3604         enddo
3605         call int_from_cart(.true.,.false.)
3606         call sc_loc_geom(.true.)
3607         do i=1,nres
3608           thetaref(i)=theta(i)
3609           phiref(i)=phi(i)
3610         enddo
3611         if (waga_dist.ne.0.0d0) then
3612           ii=0
3613           do i = nnt,nct-2
3614             do j=i+2,nct
3615
3616               x12=c(1,i)-c(1,j)
3617               y12=c(2,i)-c(2,j)
3618               z12=c(3,i)-c(3,j)
3619               distal=dsqrt(x12*x12+y12*y12+z12*z12)
3620 c              write (iout,*) k,i,j,distal,dist2_cut
3621
3622             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3623      &            .and. distal.le.dist2_cut ) then
3624
3625               ii=ii+1
3626               ii_in_use(ii)=1
3627               l_homo(k,ii)=.true.
3628
3629 c             write (iout,*) "k",k
3630 c             write (iout,*) "i",i," j",j," constr_homology",
3631 c    &                       constr_homology
3632               ires_homo(ii)=i
3633               jres_homo(ii)=j
3634               odl(k,ii)=distal
3635               if (read2sigma) then
3636                 sigma_odl(k,ii)=0
3637                 do ik=i,j
3638                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3639                 enddo
3640                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3641                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3642      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3643               else
3644                 if (odl(k,ii).le.dist_cut) then
3645                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3646                 else
3647 #ifdef OLDSIGMA
3648                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3649      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3650 #else
3651                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3652      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3653 #endif
3654                 endif
3655               endif
3656               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3657             else
3658               ii=ii+1
3659 c              l_homo(k,ii)=.false.
3660             endif
3661             enddo
3662           enddo
3663         lim_odl=ii
3664         endif
3665 c
3666 c     Theta, dihedral and SC retraints
3667 c
3668         if (waga_angle.gt.0.0d0) then
3669           do i = nnt+3,nct
3670             if (idomain(k,i).eq.0) then
3671 c               sigma_dih(k,i)=0.0
3672                cycle
3673             endif
3674             dih(k,i)=phiref(i)
3675             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3676      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
3677 c            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3678 c     &       " sigma_dihed",sigma_dih(k,i)
3679             if (sigma_dih(k,i).ne.0)
3680      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3681           enddo
3682           lim_dih=nct-nnt-2
3683         endif
3684
3685         if (waga_theta.gt.0.0d0) then
3686           do i = nnt+2,nct
3687              if (idomain(k,i).eq.0) then
3688 c              sigma_theta(k,i)=0.0
3689               cycle
3690              endif
3691              thetatpl(k,i)=thetaref(i)
3692              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3693      &                        rescore(k,i-2))/3.0
3694              if (sigma_theta(k,i).ne.0)
3695      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3696           enddo
3697         endif
3698
3699         if (waga_d.gt.0.0d0) then
3700           do i = nnt,nct
3701                if (itype(i).eq.10) cycle
3702                if (idomain(k,i).eq.0 ) then
3703 c                  sigma_d(k,i)=0.0
3704                   cycle
3705                endif
3706                xxtpl(k,i)=xxref(i)
3707                yytpl(k,i)=yyref(i)
3708                zztpl(k,i)=zzref(i)
3709                sigma_d(k,i)=rescore(k,i)
3710                if (sigma_d(k,i).ne.0)
3711      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3712                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3713           enddo
3714         endif
3715       enddo ! l
3716       enddo ! ll
3717 c
3718 c remove distance restraints not used in any model from the list
3719 c shift data in all arrays
3720 c
3721       if (waga_dist.ne.0.0d0) then
3722         ii=0
3723         liiflag=.true.
3724         do i=nnt,nct-2
3725          do j=i+2,nct
3726           ii=ii+1
3727           if (ii_in_use(ii).eq.0.and.liiflag) then
3728             liiflag=.false.
3729             iistart=ii
3730           endif
3731           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3732      &                   .not.liiflag.and.ii.eq.lim_odl) then
3733              if (ii.eq.lim_odl) then
3734               iishift=ii-iistart+1
3735              else
3736               iishift=ii-iistart
3737              endif
3738              liiflag=.true.
3739              do ki=iistart,lim_odl-iishift
3740               ires_homo(ki)=ires_homo(ki+iishift)
3741               jres_homo(ki)=jres_homo(ki+iishift)
3742               ii_in_use(ki)=ii_in_use(ki+iishift)
3743               do k=1,constr_homology
3744                odl(k,ki)=odl(k,ki+iishift)
3745                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3746                l_homo(k,ki)=l_homo(k,ki+iishift)
3747               enddo
3748              enddo
3749              ii=ii-iishift
3750              lim_odl=lim_odl-iishift
3751           endif
3752          enddo
3753         enddo
3754       endif
3755
3756       return
3757    10 stop "Error in fragment file"
3758       end