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       write (iout,*) "SETUP_VAR ialph"
1572       do i=2,nres-1
1573         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1574           nside=nside+1
1575           ialph(i,1)=nvar+nside
1576           ialph(nside,2)=i
1577         endif
1578       enddo
1579       if (indphi.gt.0) then
1580         nvar=nphi
1581       else if (indback.gt.0) then
1582         nvar=nphi+ntheta
1583       else
1584         nvar=nvar+2*nside
1585       endif
1586       write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1587       return
1588       end
1589 c----------------------------------------------------------------------------
1590       subroutine gen_dist_constr
1591 C Generate CA distance constraints.
1592       implicit none
1593       include 'DIMENSIONS'
1594       include 'COMMON.IOUNITS'
1595       include 'COMMON.GEO'
1596       include 'COMMON.VAR'
1597       include 'COMMON.INTERACT'
1598       include 'COMMON.LOCAL'
1599       include 'COMMON.NAMES'
1600       include 'COMMON.CHAIN'
1601       include 'COMMON.FFIELD'
1602       include 'COMMON.SBRIDGE'
1603       include 'COMMON.HEADER'
1604       include 'COMMON.CONTROL'
1605       include 'COMMON.DBASE'
1606       include 'COMMON.THREAD'
1607       include 'COMMON.TIME1'
1608       integer i,j,itype_pdb(maxres)
1609       common /pizda/ itype_pdb
1610       double precision dist
1611       character*2 iden
1612 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1613 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1614 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1615 cd     & ' nsup',nsup
1616       do i=nstart_sup,nstart_sup+nsup-1
1617 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1618 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1619         do j=i+2,nstart_sup+nsup-1
1620           nhpb=nhpb+1
1621           ihpb(nhpb)=i+nstart_seq-nstart_sup
1622           jhpb(nhpb)=j+nstart_seq-nstart_sup
1623           forcon(nhpb)=weidis
1624           dhpb(nhpb)=dist(i,j)
1625         enddo
1626       enddo
1627 cd      write (iout,'(a)') 'Distance constraints:'
1628 cd      do i=nss+1,nhpb
1629 cd        ii=ihpb(i)
1630 cd        jj=jhpb(i)
1631 cd        iden='CA'
1632 cd        if (ii.gt.nres) then
1633 cd          iden='SC'
1634 cd          ii=ii-nres
1635 cd          jj=jj-nres
1636 cd        endif
1637 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1638 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1639 cd     &  dhpb(i),forcon(i)
1640 cd      enddo
1641       return
1642       end
1643 c----------------------------------------------------------------------------
1644       subroutine map_read
1645       implicit none
1646       include 'DIMENSIONS'
1647       include 'COMMON.MAP'
1648       include 'COMMON.IOUNITS'
1649       integer imap
1650       character*3 angid(4) /'THE','PHI','ALP','OME'/
1651       character*80 mapcard,ucase
1652       do imap=1,nmap
1653         read (inp,'(a)') mapcard
1654         mapcard=ucase(mapcard)
1655         if (index(mapcard,'PHI').gt.0) then
1656           kang(imap)=1
1657         else if (index(mapcard,'THE').gt.0) then
1658           kang(imap)=2
1659         else if (index(mapcard,'ALP').gt.0) then
1660           kang(imap)=3
1661         else if (index(mapcard,'OME').gt.0) then
1662           kang(imap)=4
1663         else
1664           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1665           stop 'Error - illegal variable spec in MAP card.'
1666         endif
1667         call readi (mapcard,'RES1',res1(imap),0)
1668         call readi (mapcard,'RES2',res2(imap),0)
1669         if (res1(imap).eq.0) then
1670           res1(imap)=res2(imap)
1671         else if (res2(imap).eq.0) then
1672           res2(imap)=res1(imap)
1673         endif
1674         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1675           write (iout,'(a)')
1676      &    'Error - illegal definition of variable group in MAP.'
1677           stop 'Error - illegal definition of variable group in MAP.'
1678         endif
1679         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1680         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1681         call readi(mapcard,'NSTEP',nstep(imap),0)
1682         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1683           write (iout,'(a)')
1684      &     'Illegal boundary and/or step size specification in MAP.'
1685           stop 'Illegal boundary and/or step size specification in MAP.'
1686         endif
1687       enddo ! imap
1688       return
1689       end
1690 c----------------------------------------------------------------------------
1691       subroutine csaread
1692       implicit none
1693       include 'DIMENSIONS'
1694       include 'COMMON.IOUNITS'
1695       include 'COMMON.GEO'
1696       include 'COMMON.CSA'
1697       include 'COMMON.BANK'
1698       include 'COMMON.CONTROL'
1699       character*80 ucase
1700       character*620 mcmcard
1701       call card_concat(mcmcard)
1702
1703       call readi(mcmcard,'NCONF',nconf,50)
1704       call readi(mcmcard,'NADD',nadd,0)
1705       call readi(mcmcard,'JSTART',jstart,1)
1706       call readi(mcmcard,'JEND',jend,1)
1707       call readi(mcmcard,'NSTMAX',nstmax,500000)
1708       call readi(mcmcard,'N0',n0,1)
1709       call readi(mcmcard,'N1',n1,6)
1710       call readi(mcmcard,'N2',n2,4)
1711       call readi(mcmcard,'N3',n3,0)
1712       call readi(mcmcard,'N4',n4,0)
1713       call readi(mcmcard,'N5',n5,0)
1714       call readi(mcmcard,'N6',n6,10)
1715       call readi(mcmcard,'N7',n7,0)
1716       call readi(mcmcard,'N8',n8,0)
1717       call readi(mcmcard,'N9',n9,0)
1718       call readi(mcmcard,'N14',n14,0)
1719       call readi(mcmcard,'N15',n15,0)
1720       call readi(mcmcard,'N16',n16,0)
1721       call readi(mcmcard,'N17',n17,0)
1722       call readi(mcmcard,'N18',n18,0)
1723
1724       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1725
1726       call readi(mcmcard,'NDIFF',ndiff,2)
1727       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1728       call readi(mcmcard,'IS1',is1,1)
1729       call readi(mcmcard,'IS2',is2,8)
1730       call readi(mcmcard,'NRAN0',nran0,4)
1731       call readi(mcmcard,'NRAN1',nran1,2)
1732       call readi(mcmcard,'IRR',irr,1)
1733       call readi(mcmcard,'NSEED',nseed,20)
1734       call readi(mcmcard,'NTOTAL',ntotal,10000)
1735       call reada(mcmcard,'CUT1',cut1,2.0d0)
1736       call reada(mcmcard,'CUT2',cut2,5.0d0)
1737       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1738       call readi(mcmcard,'ICMAX',icmax,3)
1739       call readi(mcmcard,'IRESTART',irestart,0)
1740 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1741       ntbankm=0
1742 c!bankt
1743       call reada(mcmcard,'DELE',dele,20.0d0)
1744       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1745       call readi(mcmcard,'IREF',iref,0)
1746       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1747       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1748       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1749       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1750       write (iout,*) "NCONF_IN",nconf_in
1751       return
1752       end
1753 c----------------------------------------------------------------------------
1754       subroutine mcmread
1755       implicit none
1756       include 'DIMENSIONS'
1757       include 'COMMON.MCM'
1758       include 'COMMON.MCE'
1759       include 'COMMON.IOUNITS'
1760       character*80 ucase
1761       character*320 mcmcard
1762       integer i
1763       call card_concat(mcmcard)
1764       call readi(mcmcard,'MAXACC',maxacc,100)
1765       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1766       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1767       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1768       call readi(mcmcard,'MAXREPM',maxrepm,200)
1769       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1770       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1771       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1772       call reada(mcmcard,'E_UP',e_up,5.0D0)
1773       call reada(mcmcard,'DELTE',delte,0.1D0)
1774       call readi(mcmcard,'NSWEEP',nsweep,5)
1775       call readi(mcmcard,'NSTEPH',nsteph,0)
1776       call readi(mcmcard,'NSTEPC',nstepc,0)
1777       call reada(mcmcard,'TMIN',tmin,298.0D0)
1778       call reada(mcmcard,'TMAX',tmax,298.0D0)
1779       call readi(mcmcard,'NWINDOW',nwindow,0)
1780       call readi(mcmcard,'PRINT_MC',print_mc,0)
1781       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1782       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1783       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1784       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1785       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1786       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1787       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1788       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1789       if (nwindow.gt.0) then
1790         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1791         do i=1,nwindow
1792           winlen(i)=winend(i)-winstart(i)+1
1793         enddo
1794       endif
1795       if (tmax.lt.tmin) tmax=tmin
1796       if (tmax.eq.tmin) then
1797         nstepc=0
1798         nsteph=0
1799       endif
1800       if (nstepc.gt.0 .and. nsteph.gt.0) then
1801         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1802         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1803       endif
1804 C Probabilities of different move types
1805       sumpro_type(0)=0.0D0
1806       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1807       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1808       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1809       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1810       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1811       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1812       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1813       do i=1,MaxMoveType
1814         print *,'i',i,' sumprotype',sumpro_type(i)
1815         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1816         print *,'i',i,' sumprotype',sumpro_type(i)
1817       enddo
1818       return
1819       end
1820 c----------------------------------------------------------------------------
1821       subroutine read_minim
1822       implicit none
1823       include 'DIMENSIONS'
1824       include 'COMMON.MINIM'
1825       include 'COMMON.IOUNITS'
1826       character*80 ucase
1827       character*320 minimcard
1828       call card_concat(minimcard)
1829       call readi(minimcard,'MAXMIN',maxmin,2000)
1830       call readi(minimcard,'MAXFUN',maxfun,5000)
1831       call readi(minimcard,'MINMIN',minmin,maxmin)
1832       call readi(minimcard,'MINFUN',minfun,maxmin)
1833       call reada(minimcard,'TOLF',tolf,1.0D-2)
1834       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1835       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1836       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1837       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1838       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1839      &         'Options in energy minimization:'
1840       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1841      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1842      & 'MinMin:',MinMin,' MinFun:',MinFun,
1843      & ' TolF:',TolF,' RTolF:',RTolF
1844       return
1845       end
1846 c----------------------------------------------------------------------------
1847       subroutine read_angles(kanal,*)
1848       implicit none
1849       include 'DIMENSIONS'
1850       include 'COMMON.GEO'
1851       include 'COMMON.VAR'
1852       include 'COMMON.CHAIN'
1853       include 'COMMON.IOUNITS'
1854       include 'COMMON.CONTROL'
1855       integer i,kanal
1856 c Read angles from input
1857 c
1858        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1859        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1860        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1861        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1862
1863        do i=1,nres
1864 c 9/7/01 avoid 180 deg valence angle
1865         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1866 c
1867         theta(i)=deg2rad*theta(i)
1868         phi(i)=deg2rad*phi(i)
1869         alph(i)=deg2rad*alph(i)
1870         omeg(i)=deg2rad*omeg(i)
1871        enddo
1872       return
1873    10 return1
1874       end
1875 c----------------------------------------------------------------------------
1876       subroutine reada(rekord,lancuch,wartosc,default)
1877       implicit none
1878       character*(*) rekord,lancuch
1879       double precision wartosc,default
1880       integer ilen,iread
1881       external ilen
1882       iread=index(rekord,lancuch)
1883       if (iread.eq.0) then
1884         wartosc=default
1885         return
1886       endif
1887       iread=iread+ilen(lancuch)+1
1888       read (rekord(iread:),*,err=10,end=10) wartosc
1889       return
1890   10  wartosc=default
1891       return
1892       end
1893 c----------------------------------------------------------------------------
1894       subroutine readi(rekord,lancuch,wartosc,default)
1895       implicit none
1896       character*(*) rekord,lancuch
1897       integer wartosc,default
1898       integer ilen,iread
1899       external ilen
1900       iread=index(rekord,lancuch)
1901       if (iread.eq.0) then
1902         wartosc=default
1903         return
1904       endif
1905       iread=iread+ilen(lancuch)+1
1906       read (rekord(iread:),*,err=10,end=10) wartosc
1907       return
1908   10  wartosc=default
1909       return
1910       end
1911 c----------------------------------------------------------------------------
1912       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1913       implicit none
1914       integer dim,i
1915       integer tablica(dim),default
1916       character*(*) rekord,lancuch
1917       character*80 aux
1918       integer ilen,iread
1919       external ilen
1920       do i=1,dim
1921         tablica(i)=default
1922       enddo
1923       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1924       if (iread.eq.0) return
1925       iread=iread+ilen(lancuch)+1
1926       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1927    10 return
1928       end
1929 c----------------------------------------------------------------------------
1930       subroutine multreada(rekord,lancuch,tablica,dim,default)
1931       implicit none
1932       integer dim,i
1933       double precision tablica(dim),default
1934       character*(*) rekord,lancuch
1935       character*80 aux
1936       integer ilen,iread
1937       external ilen
1938       do i=1,dim
1939         tablica(i)=default
1940       enddo
1941       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1942       if (iread.eq.0) return
1943       iread=iread+ilen(lancuch)+1
1944       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1945    10 return
1946       end
1947 c----------------------------------------------------------------------------
1948       subroutine openunits
1949       implicit none
1950       include 'DIMENSIONS'
1951 #ifdef MPI
1952       include 'mpif.h'
1953       integer ierror
1954       character*16 form,nodename
1955       integer nodelen
1956 #endif
1957       include 'COMMON.SETUP'
1958       include 'COMMON.IOUNITS'
1959       include 'COMMON.MD'
1960       include 'COMMON.CONTROL'
1961       integer lenpre,lenpot,ilen,lentmp,npos
1962       external ilen
1963       character*3 out1file_text,ucase
1964       character*3 ll
1965       external ucase
1966       tmpdir=""
1967       out1file_text="YES"
1968 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1969       call getenv_loc("PREFIX",prefix)
1970       pref_orig = prefix
1971       call getenv_loc("POT",pot)
1972       call getenv_loc("DIRTMP",tmpdir)
1973       call getenv_loc("CURDIR",curdir)
1974       call getenv_loc("OUT1FILE",out1file_text)
1975 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1976       out1file_text=ucase(out1file_text)
1977       if (out1file_text(1:1).eq."Y") then
1978         out1file=.true.
1979       else
1980         out1file=fg_rank.gt.0
1981       endif
1982       lenpre=ilen(prefix)
1983       lenpot=ilen(pot)
1984       lentmp=ilen(tmpdir)
1985       if (lentmp.gt.0) then
1986           write (*,'(80(1h!))')
1987           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1988           write (*,'(80(1h!))')
1989           write (*,*)"All output files will be on node /tmp directory."
1990 #ifdef MPI
1991         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1992         if (me.eq.king) then
1993           write (*,*) "The master node is ",nodename
1994         else if (fg_rank.eq.0) then
1995           write (*,*) "I am the CG slave node ",nodename
1996         else
1997           write (*,*) "I am the FG slave node ",nodename
1998         endif
1999 #endif
2000         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2001         lenpre = lentmp+lenpre+1
2002       endif
2003       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2004 C Get the names and open the input files
2005 #if defined(WINIFL) || defined(WINPGI)
2006       open(1,file=pref_orig(:ilen(pref_orig))//
2007      &  '.inp',status='old',readonly,shared)
2008        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2009 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2010 C Get parameter filenames and open the parameter files.
2011       call getenv_loc('BONDPAR',bondname)
2012       open (ibond,file=bondname,status='old',readonly,shared)
2013       call getenv_loc('THETPAR',thetname)
2014       open (ithep,file=thetname,status='old',readonly,shared)
2015       call getenv_loc('ROTPAR',rotname)
2016       open (irotam,file=rotname,status='old',readonly,shared)
2017       call getenv_loc('TORPAR',torname)
2018       open (itorp,file=torname,status='old',readonly,shared)
2019       call getenv_loc('TORDPAR',tordname)
2020       open (itordp,file=tordname,status='old',readonly,shared)
2021       call getenv_loc('FOURIER',fouriername)
2022       open (ifourier,file=fouriername,status='old',readonly,shared)
2023       call getenv_loc('ELEPAR',elename)
2024       open (ielep,file=elename,status='old',readonly,shared)
2025       call getenv_loc('SIDEPAR',sidename)
2026       open (isidep,file=sidename,status='old',readonly,shared)
2027       call getenv_loc('LIPTRANPAR',liptranname)
2028       open (iliptranpar,file=liptranname,status='old',readonly,shared)
2029 #elif (defined CRAY) || (defined AIX)
2030       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2031      &  action='read')
2032 c      print *,"Processor",myrank," opened file 1"
2033       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2034 c      print *,"Processor",myrank," opened file 9"
2035 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2036 C Get parameter filenames and open the parameter files.
2037       call getenv_loc('BONDPAR',bondname)
2038       open (ibond,file=bondname,status='old',action='read')
2039 c      print *,"Processor",myrank," opened file IBOND"
2040       call getenv_loc('THETPAR',thetname)
2041       open (ithep,file=thetname,status='old',action='read')
2042 #ifndef CRYST_THETA
2043       call getenv_loc('THETPARPDB',thetname_pdb)
2044       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2045 #endif
2046 c      print *,"Processor",myrank," opened file ITHEP"
2047       call getenv_loc('ROTPAR',rotname)
2048       open (irotam,file=rotname,status='old',action='read')
2049 #ifndef CRYST_SC
2050       call getenv_loc('ROTPARPDB',rotname_pdb)
2051       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2052 #endif
2053 c      print *,"Processor",myrank," opened file IROTAM"
2054       call getenv_loc('TORPAR',torname)
2055       open (itorp,file=torname,status='old',action='read')
2056 c      print *,"Processor",myrank," opened file ITORP"
2057       call getenv_loc('TORDPAR',tordname)
2058       open (itordp,file=tordname,status='old',action='read')
2059 c      print *,"Processor",myrank," opened file ITORDP"
2060       call getenv_loc('SCCORPAR',sccorname)
2061       open (isccor,file=sccorname,status='old',action='read')
2062 c      print *,"Processor",myrank," opened file ISCCOR"
2063       call getenv_loc('FOURIER',fouriername)
2064       open (ifourier,file=fouriername,status='old',action='read')
2065 c      print *,"Processor",myrank," opened file IFOURIER"
2066       call getenv_loc('ELEPAR',elename)
2067       open (ielep,file=elename,status='old',action='read')
2068 c      print *,"Processor",myrank," opened file IELEP"
2069       call getenv_loc('SIDEPAR',sidename)
2070       open (isidep,file=sidename,status='old',action='read')
2071       call getenv_loc('LIPTRANPAR',liptranname)
2072       open (iliptranpar,file=liptranname,status='old',action='read')
2073 c      print *,"Processor",myrank," opened file ISIDEP"
2074 c      print *,"Processor",myrank," opened parameter files"
2075 #elif (defined G77)
2076       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2077       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2078 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2079 C Get parameter filenames and open the parameter files.
2080       call getenv_loc('BONDPAR',bondname)
2081       open (ibond,file=bondname,status='old')
2082       call getenv_loc('THETPAR',thetname)
2083       open (ithep,file=thetname,status='old')
2084 #ifndef CRYST_THETA
2085       call getenv_loc('THETPARPDB',thetname_pdb)
2086       open (ithep_pdb,file=thetname_pdb,status='old')
2087 #endif
2088       call getenv_loc('ROTPAR',rotname)
2089       open (irotam,file=rotname,status='old')
2090 #ifndef CRYST_SC
2091       call getenv_loc('ROTPARPDB',rotname_pdb)
2092       open (irotam_pdb,file=rotname_pdb,status='old')
2093 #endif
2094       call getenv_loc('TORPAR',torname)
2095       open (itorp,file=torname,status='old')
2096       call getenv_loc('TORDPAR',tordname)
2097       open (itordp,file=tordname,status='old')
2098       call getenv_loc('SCCORPAR',sccorname)
2099       open (isccor,file=sccorname,status='old')
2100       call getenv_loc('FOURIER',fouriername)
2101       open (ifourier,file=fouriername,status='old')
2102       call getenv_loc('ELEPAR',elename)
2103       open (ielep,file=elename,status='old')
2104       call getenv_loc('SIDEPAR',sidename)
2105       open (isidep,file=sidename,status='old')
2106       call getenv_loc('LIPTRANPAR',liptranname)
2107       open (iliptranpar,file=liptranname,status='old')
2108 #else
2109       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2110      &  readonly)
2111        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2112 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2113 C Get parameter filenames and open the parameter files.
2114       call getenv_loc('BONDPAR',bondname)
2115       open (ibond,file=bondname,status='old',readonly)
2116       call getenv_loc('THETPAR',thetname)
2117       open (ithep,file=thetname,status='old',readonly)
2118       call getenv_loc('ROTPAR',rotname)
2119       open (irotam,file=rotname,status='old',readonly)
2120       call getenv_loc('TORPAR',torname)
2121       open (itorp,file=torname,status='old',readonly)
2122       call getenv_loc('TORDPAR',tordname)
2123       open (itordp,file=tordname,status='old',readonly)
2124       call getenv_loc('SCCORPAR',sccorname)
2125       open (isccor,file=sccorname,status='old',readonly)
2126 #ifndef CRYST_THETA
2127       call getenv_loc('THETPARPDB',thetname_pdb)
2128       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2129 #endif
2130       call getenv_loc('FOURIER',fouriername)
2131       open (ifourier,file=fouriername,status='old',readonly)
2132       call getenv_loc('ELEPAR',elename)
2133       open (ielep,file=elename,status='old',readonly)
2134       call getenv_loc('SIDEPAR',sidename)
2135       open (isidep,file=sidename,status='old',readonly)
2136       call getenv_loc('LIPTRANPAR',liptranname)
2137       open (iliptranpar,file=liptranname,status='old',action='read')
2138 #ifndef CRYST_SC
2139       call getenv_loc('ROTPARPDB',rotname_pdb)
2140       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2141 #endif
2142 #endif
2143 #ifdef TUBE
2144       call getenv_loc('TUBEPAR',tubename)
2145 #if defined(WINIFL) || defined(WINPGI)
2146       open (itube,file=tubename,status='old',readonly,shared)
2147 #elif (defined CRAY)  || (defined AIX)
2148       open (itube,file=tubename,status='old',action='read')
2149 #elif (defined G77)
2150       open (itube,file=tubename,status='old')
2151 #else
2152       open (itube,file=tubename,status='old',readonly)
2153 #endif
2154 #endif
2155 #ifndef OLDSCP
2156 C
2157 C 8/9/01 In the newest version SCp interaction constants are read from a file
2158 C Use -DOLDSCP to use hard-coded constants instead.
2159 C
2160       call getenv_loc('SCPPAR',scpname)
2161 #if defined(WINIFL) || defined(WINPGI)
2162       open (iscpp,file=scpname,status='old',readonly,shared)
2163 #elif (defined CRAY)  || (defined AIX)
2164       open (iscpp,file=scpname,status='old',action='read')
2165 #elif (defined G77)
2166       open (iscpp,file=scpname,status='old')
2167 #else
2168       open (iscpp,file=scpname,status='old',readonly)
2169 #endif
2170 #endif
2171       call getenv_loc('PATTERN',patname)
2172 #if defined(WINIFL) || defined(WINPGI)
2173       open (icbase,file=patname,status='old',readonly,shared)
2174 #elif (defined CRAY)  || (defined AIX)
2175       open (icbase,file=patname,status='old',action='read')
2176 #elif (defined G77)
2177       open (icbase,file=patname,status='old')
2178 #else
2179       open (icbase,file=patname,status='old',readonly)
2180 #endif
2181 #ifdef MPI
2182 C Open output file only for CG processes
2183 c      print *,"Processor",myrank," fg_rank",fg_rank
2184       if (fg_rank.eq.0) then
2185
2186       if (nodes.eq.1) then
2187         npos=3
2188       else
2189         npos = dlog10(dfloat(nodes-1))+1
2190       endif
2191       if (npos.lt.3) npos=3
2192       write (liczba,'(i1)') npos
2193       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2194      &  //')'
2195       write (liczba,form) me
2196       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2197      &  liczba(:ilen(liczba))
2198       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2199      &  //'.int'
2200       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2201      &  //'.pdb'
2202       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2203      &  liczba(:ilen(liczba))//'.mol2'
2204       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2205      &  liczba(:ilen(liczba))//'.stat'
2206       if (lentmp.gt.0)
2207      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2208      &      //liczba(:ilen(liczba))//'.stat')
2209       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2210      &  //'.rst'
2211       if(usampl) then
2212           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2213      & liczba(:ilen(liczba))//'.const'
2214       endif
2215
2216       endif
2217 #else
2218       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2219       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2220       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2221       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2222       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2223       if (lentmp.gt.0)
2224      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2225      &    //'.stat')
2226       rest2name=prefix(:ilen(prefix))//'.rst'
2227       if(usampl) then
2228          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2229       endif
2230 #endif
2231 #if defined(AIX) || defined(PGI) || defined(CRAY)
2232       if (me.eq.king .or. .not. out1file) then
2233          open(iout,file=outname,status='unknown')
2234       else
2235          open(iout,file="/dev/null",status="unknown")
2236       endif
2237 c#define DEBUG
2238 #ifdef DEBUG
2239       if (fg_rank.gt.0) then
2240         write (liczba,'(i3.3)') myrank/nfgtasks
2241         write (ll,'(bz,i3.3)') fg_rank
2242         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2243      &   status='unknown')
2244       endif
2245 #endif
2246 c#undef DEBUG
2247       if(me.eq.king) then
2248        open(igeom,file=intname,status='unknown',position='append')
2249        open(ipdb,file=pdbname,status='unknown')
2250        open(imol2,file=mol2name,status='unknown')
2251        open(istat,file=statname,status='unknown',position='append')
2252       else
2253 c1out       open(iout,file=outname,status='unknown')
2254       endif
2255 #else
2256       if (me.eq.king .or. .not.out1file)
2257      &    open(iout,file=outname,status='unknown')
2258 #ifdef DEBUG
2259       if (fg_rank.gt.0) then
2260         write (liczba,'(i3.3)') myrank/nfgtasks
2261         write (ll,'(bz,i3.3)') fg_rank
2262         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2263      &   status='unknown')
2264       endif
2265 #endif
2266       if(me.eq.king) then
2267        open(igeom,file=intname,status='unknown',access='append')
2268        open(ipdb,file=pdbname,status='unknown')
2269        open(imol2,file=mol2name,status='unknown')
2270        open(istat,file=statname,status='unknown',access='append')
2271       else
2272 c1out       open(iout,file=outname,status='unknown')
2273       endif
2274 #endif
2275       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2276       csa_seed=prefix(:lenpre)//'.CSA.seed'
2277       csa_history=prefix(:lenpre)//'.CSA.history'
2278       csa_bank=prefix(:lenpre)//'.CSA.bank'
2279       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2280       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2281       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2282 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2283       csa_int=prefix(:lenpre)//'.int'
2284       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2285       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2286       csa_in=prefix(:lenpre)//'.CSA.in'
2287 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2288 C Write file names
2289       if (me.eq.king)then
2290       write (iout,'(80(1h-))')
2291       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2292       write (iout,'(80(1h-))')
2293       write (iout,*) "Input file                      : ",
2294      &  pref_orig(:ilen(pref_orig))//'.inp'
2295       write (iout,*) "Output file                     : ",
2296      &  outname(:ilen(outname))
2297       write (iout,*)
2298       write (iout,*) "Sidechain potential file        : ",
2299      &  sidename(:ilen(sidename))
2300 #ifndef OLDSCP
2301       write (iout,*) "SCp potential file              : ",
2302      &  scpname(:ilen(scpname))
2303 #endif
2304       write (iout,*) "Electrostatic potential file    : ",
2305      &  elename(:ilen(elename))
2306       write (iout,*) "Cumulant coefficient file       : ",
2307      &  fouriername(:ilen(fouriername))
2308       write (iout,*) "Torsional parameter file        : ",
2309      &  torname(:ilen(torname))
2310       write (iout,*) "Double torsional parameter file : ",
2311      &  tordname(:ilen(tordname))
2312       write (iout,*) "SCCOR parameter file : ",
2313      &  sccorname(:ilen(sccorname))
2314       write (iout,*) "Bond & inertia constant file    : ",
2315      &  bondname(:ilen(bondname))
2316       write (iout,*) "Bending-torsion parameter file  : ",
2317      &  thetname(:ilen(thetname))
2318       write (iout,*) "Rotamer parameter file          : ",
2319      &  rotname(:ilen(rotname))
2320       write (iout,*) "Thetpdb parameter file          : ",
2321      &  thetname_pdb(:ilen(thetname_pdb))
2322       write (iout,*) "Threading database              : ",
2323      &  patname(:ilen(patname))
2324       if (lentmp.ne.0)
2325      &write (iout,*)" DIRTMP                          : ",
2326      &  tmpdir(:lentmp)
2327       write (iout,'(80(1h-))')
2328       endif
2329       return
2330       end
2331 c----------------------------------------------------------------------------
2332       subroutine card_concat(card)
2333       implicit real*8 (a-h,o-z)
2334       include 'DIMENSIONS'
2335       include 'COMMON.IOUNITS'
2336       character*(*) card
2337       character*80 karta,ucase
2338       external ilen
2339       read (inp,'(a)') karta
2340       karta=ucase(karta)
2341       card=' '
2342       do while (karta(80:80).eq.'&')
2343         card=card(:ilen(card)+1)//karta(:79)
2344         read (inp,'(a)') karta
2345         karta=ucase(karta)
2346       enddo
2347       card=card(:ilen(card)+1)//karta
2348       return
2349       end
2350 c------------------------------------------------------------------------------
2351       subroutine readrst
2352       implicit none
2353       include 'DIMENSIONS'
2354       include 'COMMON.CHAIN'
2355       include 'COMMON.IOUNITS'
2356       include 'COMMON.CONTROL'
2357       include 'COMMON.MD'
2358       include 'COMMON.QRESTR'
2359 #ifdef FIVEDIAG
2360        include 'COMMON.LAGRANGE.5diag'
2361 #else
2362        include 'COMMON.LAGRANGE'
2363 #endif
2364       integer i,j
2365       open(irest2,file=rest2name,status='unknown')
2366       read(irest2,*) totT,EK,potE,totE,t_bath
2367       totTafm=totT
2368       do i=1,2*nres
2369          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2370       enddo
2371       do i=1,2*nres
2372          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2373       enddo
2374       if(usampl) then
2375              read (irest2,*) iset
2376       endif
2377       close(irest2)
2378       return
2379       end
2380 c------------------------------------------------------------------------------
2381       subroutine read_fragments
2382       implicit none
2383       include 'DIMENSIONS'
2384 #ifdef MPI
2385       include 'mpif.h'
2386 #endif
2387       include 'COMMON.SETUP'
2388       include 'COMMON.CHAIN'
2389       include 'COMMON.IOUNITS'
2390       include 'COMMON.MD'
2391       include 'COMMON.QRESTR'
2392       include 'COMMON.CONTROL'
2393       integer i
2394       read(inp,*) nset,nfrag,npair,nfrag_back
2395       loc_qlike=(nfrag_back.lt.0)
2396       nfrag_back=iabs(nfrag_back)
2397       if(me.eq.king.or..not.out1file)
2398      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2399      &  " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2400       do iset=1,nset
2401          read(inp,*) mset(iset)
2402        do i=1,nfrag
2403          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2404      &     qinfrag(i,iset)
2405          if(me.eq.king.or..not.out1file)
2406      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2407      &     ifrag(2,i,iset), qinfrag(i,iset)
2408        enddo
2409        do i=1,npair
2410         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2411      &    qinpair(i,iset)
2412         if(me.eq.king.or..not.out1file)
2413      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2414      &    ipair(2,i,iset), qinpair(i,iset)
2415        enddo
2416        if (loc_qlike) then
2417        do i=1,nfrag_back
2418         read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2419      &     wfrag_back(2,i,iset),qin_back(2,i,iset),
2420      &     wfrag_back(3,i,iset),qin_back(3,i,iset),
2421      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2422         if(me.eq.king.or..not.out1file)
2423      &   write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2424      &   wfrag_back(2,i,iset),qin_back(3,i,iset),
2425      &   wfrag_back(3,i,iset),qin_back(3,i,iset),
2426      &   ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2427        enddo
2428        else
2429        do i=1,nfrag_back
2430         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2431      &     wfrag_back(3,i,iset),
2432      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2433         if(me.eq.king.or..not.out1file)
2434      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2435      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2436        enddo
2437        endif
2438       enddo
2439       return
2440       end
2441 C---------------------------------------------------------------------------
2442       subroutine read_afminp
2443       implicit none
2444       include 'DIMENSIONS'
2445 #ifdef MPI
2446       include 'mpif.h'
2447 #endif
2448       include 'COMMON.SETUP'
2449       include 'COMMON.CONTROL'
2450       include 'COMMON.CHAIN'
2451       include 'COMMON.IOUNITS'
2452       include 'COMMON.SBRIDGE'
2453       character*320 afmcard
2454       integer i
2455 c      print *, "wchodze"
2456       call card_concat(afmcard)
2457       call readi(afmcard,"BEG",afmbeg,0)
2458       call readi(afmcard,"END",afmend,0)
2459       call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2460       call reada(afmcard,"VEL",velAFMconst,0.0d0)
2461       print *,'FORCE=' ,forceAFMconst
2462 CCCC NOW PROPERTIES FOR AFM
2463        distafminit=0.0d0
2464        do i=1,3
2465         distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2466        enddo
2467         distafminit=dsqrt(distafminit)
2468 c        print *,'initdist',distafminit
2469       return
2470       end
2471 c-------------------------------------------------------------------------------
2472       subroutine read_saxs_constr
2473       implicit none
2474       include 'DIMENSIONS'
2475 #ifdef MPI
2476       include 'mpif.h'
2477 #endif
2478       include 'COMMON.SETUP'
2479       include 'COMMON.CONTROL'
2480       include 'COMMON.SAXS'
2481       include 'COMMON.CHAIN'
2482       include 'COMMON.IOUNITS'
2483       include 'COMMON.SBRIDGE'
2484       double precision cm(3),cnorm
2485       integer i,j
2486 c      read(inp,*) nsaxs
2487       write (iout,*) "Calling read_saxs nsaxs",nsaxs
2488       call flush(iout)
2489       if (saxs_mode.eq.0) then
2490 c SAXS distance distribution
2491       do i=1,nsaxs
2492         read(inp,*) distsaxs(i),Psaxs(i)
2493       enddo
2494       Cnorm = 0.0d0
2495       do i=1,nsaxs
2496         Cnorm = Cnorm + Psaxs(i)
2497       enddo
2498       write (iout,*) "Cnorm",Cnorm
2499       do i=1,nsaxs
2500         Psaxs(i)=Psaxs(i)/Cnorm
2501       enddo
2502       write (iout,*) "Normalized distance distribution from SAXS"
2503       do i=1,nsaxs
2504         write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2505       enddo
2506       Wsaxs0=0.0d0
2507       do i=1,nsaxs
2508         Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2509       enddo
2510       write (iout,*) "Wsaxs0",Wsaxs0
2511       else
2512 c SAXS "spheres".
2513       do i=1,nsaxs
2514         read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2515       enddo
2516       do j=1,3
2517         cm(j)=0.0d0
2518       enddo
2519       do i=1,nsaxs
2520         do j=1,3
2521           cm(j)=cm(j)+Csaxs(j,i)
2522         enddo
2523       enddo
2524       do j=1,3
2525         cm(j)=cm(j)/nsaxs
2526       enddo
2527       do i=1,nsaxs
2528         do j=1,3
2529           Csaxs(j,i)=Csaxs(j,i)-cm(j)
2530         enddo
2531       enddo
2532       write (iout,*) "SAXS sphere coordinates"
2533       do i=1,nsaxs
2534         write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2535       enddo
2536       endif
2537       return
2538       end
2539
2540 c-------------------------------------------------------------------------------
2541       subroutine read_dist_constr
2542       implicit none
2543       include 'DIMENSIONS'
2544 #ifdef MPI
2545       include 'mpif.h'
2546 #endif
2547       include 'COMMON.SETUP'
2548       include 'COMMON.CONTROL'
2549       include 'COMMON.CHAIN'
2550       include 'COMMON.IOUNITS'
2551       include 'COMMON.SBRIDGE'
2552       include 'COMMON.INTERACT'
2553       integer i,j,k,ii,jj,itemp,link_type,iiend,jjend,kk
2554       integer nfrag_,npair_,ndist_,ifrag_(2,100),ipair_(2,1000)
2555       double precision wfrag_(100),wpair_(1000)
2556       double precision ddjk,dist,dist_cut,fordepthmax
2557       character*5000 controlcard
2558       logical normalize,next
2559       integer restr_type
2560       double precision scal_bfac
2561       double precision xlink(4,0:4) /
2562 c           a          b       c     sigma
2563      &   0.0d0,0.0d0,0.0d0,0.0d0,                             ! default, no xlink potential
2564      &   0.00305218d0,9.46638d0,4.68901d0,4.74347d0,          ! ZL
2565      &   0.00214928d0,12.7517d0,0.00375009d0,6.13477d0,       ! ADH
2566      &   0.00184547d0,11.2678d0,0.00140292d0,7.00868d0,       ! PDH
2567      &   0.000161786d0,6.29273d0,4.40993d0,7.13956d0    /     ! DSS
2568 c      print *, "WCHODZE"
2569       write (iout,*) "Calling read_dist_constr"
2570 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2571 c      call flush(iout)
2572       restr_on_coord=.false.
2573       next=.true.
2574
2575       npeak=0
2576       ipeak=0
2577       nhpb_peak=0
2578
2579       DO WHILE (next)
2580
2581       call card_concat(controlcard)
2582       next = index(controlcard,"NEXT").gt.0
2583       call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2584       write (iout,*) "restr_type",restr_type
2585       call readi(controlcard,"NFRAG",nfrag_,0)
2586       call readi(controlcard,"NFRAG",nfrag_,0)
2587       call readi(controlcard,"NPAIR",npair_,0)
2588       call readi(controlcard,"NDIST",ndist_,0)
2589       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2590       call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2591       if (restr_type.eq.10)
2592      &  call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2593       if (restr_type.eq.12)
2594      &  call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2595       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2596       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2597       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2598       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2599       normalize = index(controlcard,"NORMALIZE").gt.0
2600       write (iout,*) "WBOLTZD",wboltzd
2601       write (iout,*) "SCAL_PEAK",scal_peak
2602       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2603       write (iout,*) "IFRAG"
2604       do i=1,nfrag_
2605         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2606       enddo
2607       write (iout,*) "IPAIR"
2608       do i=1,npair_
2609         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2610       enddo
2611       if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2612      & write (iout,*)
2613      &   "Distance restraints as generated from reference structure"
2614       do i=1,nfrag_
2615         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2616         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2617      &    ifrag_(2,i)=nstart_sup+nsup-1
2618 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2619 c        call flush(iout)
2620         if (wfrag_(i).eq.0.0d0) cycle
2621         do j=ifrag_(1,i),ifrag_(2,i)-1
2622           do k=j+1,ifrag_(2,i)
2623 c            write (iout,*) "j",j," k",k
2624             ddjk=dist(j,k)
2625             if (restr_type.eq.1) then
2626               nhpb=nhpb+1
2627               irestr_type(nhpb)=1
2628               ihpb(nhpb)=j
2629               jhpb(nhpb)=k
2630               dhpb(nhpb)=ddjk
2631               forcon(nhpb)=wfrag_(i)
2632             else if (constr_dist.eq.2) then
2633               if (ddjk.le.dist_cut) then
2634                 nhpb=nhpb+1
2635                 irestr_type(nhpb)=1
2636                 ihpb(nhpb)=j
2637                 jhpb(nhpb)=k
2638                 dhpb(nhpb)=ddjk
2639                 forcon(nhpb)=wfrag_(i)
2640               endif
2641             else if (restr_type.eq.3) then
2642               nhpb=nhpb+1
2643               irestr_type(nhpb)=1
2644               ihpb(nhpb)=j
2645               jhpb(nhpb)=k
2646               dhpb(nhpb)=ddjk
2647               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2648             endif
2649 #ifdef MPI
2650             if (.not.out1file .or. me.eq.king)
2651      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2652      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2653 #else
2654             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2655      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2656 #endif
2657           enddo
2658         enddo
2659       enddo
2660       do i=1,npair_
2661         if (wpair_(i).eq.0.0d0) cycle
2662         ii = ipair_(1,i)
2663         jj = ipair_(2,i)
2664         if (ii.gt.jj) then
2665           itemp=ii
2666           ii=jj
2667           jj=itemp
2668         endif
2669         do j=ifrag_(1,ii),ifrag_(2,ii)
2670           do k=ifrag_(1,jj),ifrag_(2,jj)
2671             ddjk=dist(j,k)
2672             if (restr_type.eq.1) then
2673               nhpb=nhpb+1
2674               irestr_type(nhpb)=1
2675               ihpb(nhpb)=j
2676               jhpb(nhpb)=k
2677               dhpb(nhpb)=ddjk
2678               forcon(nhpb)=wpair_(i)
2679             else if (constr_dist.eq.2) then
2680               if (ddjk.le.dist_cut) then
2681                 nhpb=nhpb+1
2682                 irestr_type(nhpb)=1
2683                 ihpb(nhpb)=j
2684                 jhpb(nhpb)=k
2685                 dhpb(nhpb)=ddjk
2686                 forcon(nhpb)=wpair_(i)
2687               endif
2688             else if (restr_type.eq.3) then
2689               nhpb=nhpb+1
2690               irestr_type(nhpb)=1
2691               ihpb(nhpb)=j
2692               jhpb(nhpb)=k
2693               dhpb(nhpb)=ddjk
2694               forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2695             endif
2696 #ifdef MPI
2697             if (.not.out1file .or. me.eq.king)
2698      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2699      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2700 #else
2701             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2702      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2703 #endif
2704           enddo
2705         enddo
2706       enddo
2707
2708 c      print *,ndist_
2709       write (iout,*) "Distance restraints as read from input"
2710       do i=1,ndist_
2711         if (restr_type.eq.12) then
2712           read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2713      &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2714      &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2715      &    fordepth_peak(nhpb_peak+1),npeak
2716 c          write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2717 c     &    dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2718 c     &    ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2719 c     &    fordepth_peak(nhpb_peak+1),npeak
2720           if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2721      &      fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2722           nhpb_peak=nhpb_peak+1
2723           irestr_type_peak(nhpb_peak)=12
2724           if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2725           ipeak(2,npeak)=i
2726 #ifdef MPI
2727           if (.not.out1file .or. me.eq.king)
2728      &    write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2729      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2730      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2731      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2732      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2733 #else
2734           write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2735      &     nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2736      &     ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2737      &     dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2738      &     fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2739 #endif
2740           if (ibecarb_peak(nhpb_peak).eq.3) then
2741             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2742           else if (ibecarb_peak(nhpb_peak).eq.2) then
2743             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2744           else if (ibecarb_peak(nhpb_peak).eq.1) then
2745             ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2746             jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2747           endif
2748         else if (restr_type.eq.11) then
2749           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2750      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2751 c        fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2752           if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2753           nhpb=nhpb+1
2754           irestr_type(nhpb)=11
2755 #ifdef MPI
2756           if (.not.out1file .or. me.eq.king)
2757      &    write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2758      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2759      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2760 #else
2761           write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2762      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2763      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2764 #endif
2765 c          if (ibecarb(nhpb).gt.0) then
2766 c            ihpb(nhpb)=ihpb(nhpb)+nres
2767 c            jhpb(nhpb)=jhpb(nhpb)+nres
2768 c          endif
2769          if (ibecarb(nhpb).eq.3) then
2770             ihpb(nhpb)=ihpb(nhpb)+nres
2771           else if (ibecarb(nhpb).eq.2) then
2772             ihpb(nhpb)=ihpb(nhpb)+nres
2773           else if (ibecarb(nhpb).eq.1) then
2774             ihpb(nhpb)=ihpb(nhpb)+nres
2775             jhpb(nhpb)=jhpb(nhpb)+nres
2776           endif
2777         else if (restr_type.eq.10) then
2778 c Cross-lonk Markov-like potential
2779           call card_concat(controlcard)
2780           call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2781           call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2782           ibecarb(nhpb+1)=0
2783           if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2784           if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2785           if (index(controlcard,"ZL").gt.0) then
2786             link_type=1
2787           else if (index(controlcard,"ADH").gt.0) then
2788             link_type=2
2789           else if (index(controlcard,"PDH").gt.0) then
2790             link_type=3
2791           else if (index(controlcard,"DSS").gt.0) then
2792             link_type=4
2793           else
2794             link_type=0
2795           endif
2796           call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2797      &       xlink(1,link_type))
2798           call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2799      &       xlink(2,link_type))
2800           call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2801      &       xlink(3,link_type))
2802           call reada(controlcard,"SIGMA",forcon(nhpb+1),
2803      &       xlink(4,link_type))
2804           call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2805 c          read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2806 c     &      dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2807           if (forcon(nhpb+1).le.0.0d0 .or.
2808      &       (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2809           nhpb=nhpb+1
2810           irestr_type(nhpb)=10
2811           if (ibecarb(nhpb).eq.3) then
2812             jhpb(nhpb)=jhpb(nhpb)+nres
2813           else if (ibecarb(nhpb).eq.2) then
2814             ihpb(nhpb)=ihpb(nhpb)+nres
2815           else if (ibecarb(nhpb).eq.1) then
2816             ihpb(nhpb)=ihpb(nhpb)+nres
2817             jhpb(nhpb)=jhpb(nhpb)+nres
2818           endif
2819 #ifdef MPI
2820           if (.not.out1file .or. me.eq.king)
2821      &    write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2822      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2823      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2824      &     irestr_type(nhpb)
2825 #else
2826           write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2827      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2828      &     dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2829      &     irestr_type(nhpb)
2830 #endif
2831         else
2832 C        print *,"in else"
2833           read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2834      &     dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2835           if (forcon(nhpb+1).gt.0.0d0) then
2836           nhpb=nhpb+1
2837           if (dhpb1(nhpb).eq.0.0d0) then
2838             irestr_type(nhpb)=1
2839           else
2840             irestr_type(nhpb)=2
2841           endif
2842           if (ibecarb(nhpb).eq.3) then
2843             jhpb(nhpb)=jhpb(nhpb)+nres
2844           else if (ibecarb(nhpb).eq.2) then
2845             ihpb(nhpb)=ihpb(nhpb)+nres
2846           else if (ibecarb(nhpb).eq.1) then
2847             ihpb(nhpb)=ihpb(nhpb)+nres
2848             jhpb(nhpb)=jhpb(nhpb)+nres
2849           endif
2850           if (dhpb(nhpb).eq.0.0d0)
2851      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2852         endif
2853 #ifdef MPI
2854           if (.not.out1file .or. me.eq.king)
2855      &    write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2856      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2857 #else
2858           write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2859      &     nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2860 #endif
2861         endif
2862 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2863 C        if (forcon(nhpb+1).gt.0.0d0) then
2864 C          nhpb=nhpb+1
2865 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2866       enddo
2867
2868       if (restr_type.eq.4) then
2869         write (iout,*) "The BFAC array"
2870         do i=nnt,nct
2871           write (iout,'(i5,f10.5)') i,bfac(i)
2872         enddo
2873         do i=nnt,nct
2874           if (itype(i).eq.ntyp1) cycle
2875           do j=nnt,i-1
2876             if (itype(j).eq.ntyp1) cycle
2877             if (itype(i).eq.10) then
2878               iiend=0
2879             else
2880               iiend=1
2881             endif
2882             if (itype(j).eq.10) then
2883               jjend=0
2884             else
2885               jjend=1
2886             endif
2887             kk=0
2888             do ii=0,iiend
2889             do jj=0,jjend
2890             nhpb=nhpb+1
2891             irestr_type(nhpb)=1
2892             forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2893             irestr_type(nhpb)=1
2894             ibecarb(nhpb)=kk
2895             if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2896             ihpb(nhpb)=i+nres*ii
2897             jhpb(nhpb)=j+nres*jj
2898             dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2899 #ifdef MPI
2900             if (.not.out1file .or. me.eq.king) then
2901             write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2902      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2903      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2904      &       irestr_type(nhpb)
2905             endif
2906 #else
2907             write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2908      &       nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2909      &       dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2910      &       irestr_type(nhpb)
2911 #endif
2912             kk=kk+1
2913           enddo
2914           enddo
2915           enddo
2916         enddo
2917       endif
2918
2919       if (restr_type.eq.5) then
2920         restr_on_coord=.true.
2921         do i=nnt,nct
2922           if (itype(i).eq.ntyp1) cycle
2923           bfac(i)=(scal_bfac/bfac(i))**2
2924         enddo
2925       endif
2926
2927       ENDDO ! next
2928
2929       fordepthmax=0.0d0
2930       if (normalize) then
2931         do i=nss+1,nhpb
2932           if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
2933      &      fordepthmax=fordepth(i)
2934         enddo
2935         do i=nss+1,nhpb
2936           if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
2937         enddo
2938       endif
2939       return
2940       end
2941 c-------------------------------------------------------------------------------
2942       subroutine read_constr_homology
2943       implicit none
2944       include 'DIMENSIONS'
2945 #ifdef MPI
2946       include 'mpif.h'
2947 #endif
2948       include 'COMMON.SETUP'
2949       include 'COMMON.CONTROL'
2950       include 'COMMON.HOMOLOGY'
2951       include 'COMMON.CHAIN'
2952       include 'COMMON.IOUNITS'
2953       include 'COMMON.MD'
2954       include 'COMMON.QRESTR'
2955       include 'COMMON.GEO'
2956       include 'COMMON.INTERACT'
2957       include 'COMMON.NAMES'
2958 c
2959 c For new homol impl
2960 c
2961       include 'COMMON.VAR'
2962 c
2963
2964 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2965 c    &                 dist_cut
2966 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2967 c    &    sigma_odl_temp(maxres,maxres,max_template)
2968       character*2 kic2
2969       character*24 model_ki_dist, model_ki_angle
2970       character*500 controlcard
2971       integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
2972      & ik,iistart,iishift
2973       integer ilen
2974       external ilen
2975       logical liiflag
2976 c
2977 c     FP - Nov. 2014 Temporary specifications for new vars
2978 c
2979       double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
2980      &                 rescore3_tmp
2981       double precision, dimension (max_template,maxres) :: rescore
2982       double precision, dimension (max_template,maxres) :: rescore2
2983       double precision, dimension (max_template,maxres) :: rescore3
2984       double precision distal
2985       character*24 pdbfile,tpl_k_rescore
2986 c -----------------------------------------------------------------
2987 c Reading multiple PDB ref structures and calculation of retraints
2988 c not using pre-computed ones stored in files model_ki_{dist,angle}
2989 c FP (Nov., 2014)
2990 c -----------------------------------------------------------------
2991 c
2992 c
2993 c Alternative: reading from input
2994       call card_concat(controlcard)
2995       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2996       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2997       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2998       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2999       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
3000       call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
3001       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3002       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3003       start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3004       if(.not.read2sigma.and.start_from_model) then
3005           if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3006      &      write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3007           start_from_model=.false.
3008       endif
3009       if(start_from_model .and. (me.eq.king .or. .not. out1file))
3010      &    write(iout,*) 'START_FROM_MODELS is ON'
3011       if(start_from_model .and. rest) then
3012         if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3013          write(iout,*) 'START_FROM_MODELS is OFF'
3014          write(iout,*) 'remove restart keyword from input'
3015         endif
3016       endif
3017       if (homol_nset.gt.1)then
3018          call card_concat(controlcard)
3019          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3020          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3021           write(iout,*) "iset homology_weight "
3022           do i=1,homol_nset
3023            write(iout,*) i,waga_homology(i)
3024           enddo
3025          endif
3026          iset=mod(kolor,homol_nset)+1
3027       else
3028        iset=1
3029        waga_homology(1)=1.0
3030       endif
3031
3032 cd      write (iout,*) "nnt",nnt," nct",nct
3033 cd      call flush(iout)
3034
3035
3036       lim_odl=0
3037       lim_dih=0
3038 c
3039 c      write(iout,*) 'nnt=',nnt,'nct=',nct
3040 c
3041       do i = nnt,nct
3042         do k=1,constr_homology
3043          idomain(k,i)=0
3044         enddo
3045       enddo
3046
3047       ii=0
3048       do i = nnt,nct-2
3049         do j=i+2,nct
3050         ii=ii+1
3051         ii_in_use(ii)=0
3052         enddo
3053       enddo
3054
3055       if (read_homol_frag) then
3056         call read_klapaucjusz
3057       else
3058
3059       do k=1,constr_homology
3060
3061         read(inp,'(a)') pdbfile
3062         if(me.eq.king .or. .not. out1file)
3063      &    write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3064      &  pdbfile(:ilen(pdbfile))
3065         open(ipdbin,file=pdbfile,status='old',err=33)
3066         goto 34
3067   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
3068      &  pdbfile(:ilen(pdbfile))
3069         stop
3070   34    continue
3071 c        print *,'Begin reading pdb data'
3072 c
3073 c Files containing res sim or local scores (former containing sigmas)
3074 c
3075
3076         write(kic2,'(bz,i2.2)') k
3077
3078         tpl_k_rescore="template"//kic2//".sco"
3079
3080         unres_pdb=.false.
3081         if (read2sigma) then
3082           call readpdb_template(k)
3083         else
3084           call readpdb
3085         endif
3086 c
3087 c     Distance restraints
3088 c
3089 c          ... --> odl(k,ii)
3090 C Copy the coordinates from reference coordinates (?)
3091         do i=1,2*nres
3092           do j=1,3
3093             c(j,i)=cref(j,i)
3094 c           write (iout,*) "c(",j,i,") =",c(j,i)
3095           enddo
3096         enddo
3097 c
3098 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3099 c
3100 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3101           open (ientin,file=tpl_k_rescore,status='old')
3102           if (nnt.gt.1) rescore(k,1)=0.0d0
3103           do irec=nnt,nct ! loop for reading res sim
3104             if (read2sigma) then
3105              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3106      &                                rescore3_tmp,idomain_tmp
3107              i_tmp=i_tmp+nnt-1
3108              idomain(k,i_tmp)=idomain_tmp
3109              rescore(k,i_tmp)=rescore_tmp
3110              rescore2(k,i_tmp)=rescore2_tmp
3111              rescore3(k,i_tmp)=rescore3_tmp
3112              if (.not. out1file .or. me.eq.king)
3113      &        write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3114      &                      i_tmp,rescore2_tmp,rescore_tmp,
3115      &                                rescore3_tmp,idomain_tmp
3116             else
3117              idomain(k,irec)=1
3118              read (ientin,*,end=1401) rescore_tmp
3119
3120 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3121              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3122 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3123             endif
3124           enddo
3125  1401   continue
3126         close (ientin)
3127         if (waga_dist.ne.0.0d0) then
3128           ii=0
3129           do i = nnt,nct-2
3130             do j=i+2,nct
3131
3132               x12=c(1,i)-c(1,j)
3133               y12=c(2,i)-c(2,j)
3134               z12=c(3,i)-c(3,j)
3135               distal=dsqrt(x12*x12+y12*y12+z12*z12)
3136 c              write (iout,*) k,i,j,distal,dist2_cut
3137
3138             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3139      &            .and. distal.le.dist2_cut ) then
3140
3141               ii=ii+1
3142               ii_in_use(ii)=1
3143               l_homo(k,ii)=.true.
3144
3145 c             write (iout,*) "k",k
3146 c             write (iout,*) "i",i," j",j," constr_homology",
3147 c    &                       constr_homology
3148               ires_homo(ii)=i
3149               jres_homo(ii)=j
3150               odl(k,ii)=distal
3151               if (read2sigma) then
3152                 sigma_odl(k,ii)=0
3153                 do ik=i,j
3154                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3155                 enddo
3156                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3157                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3158      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3159               else
3160                 if (odl(k,ii).le.dist_cut) then
3161                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3162                 else
3163 #ifdef OLDSIGMA
3164                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3165      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3166 #else
3167                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3168      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3169 #endif
3170                 endif
3171               endif
3172               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3173             else
3174               ii=ii+1
3175               l_homo(k,ii)=.false.
3176             endif
3177             enddo
3178           enddo
3179         lim_odl=ii
3180         endif
3181 c
3182 c     Theta, dihedral and SC retraints
3183 c
3184         if (waga_angle.gt.0.0d0) then
3185 c         open (ientin,file=tpl_k_sigma_dih,status='old')
3186 c         do irec=1,maxres-3 ! loop for reading sigma_dih
3187 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3188 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3189 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3190 c    &                            sigma_dih(k,i+nnt-1)
3191 c         enddo
3192 c1402   continue
3193 c         close (ientin)
3194           do i = nnt+3,nct
3195             if (idomain(k,i).eq.0) then
3196                sigma_dih(k,i)=0.0
3197                cycle
3198             endif
3199             dih(k,i)=phiref(i) ! right?
3200 c           read (ientin,*) sigma_dih(k,i) ! original variant
3201 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
3202 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3203 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
3204 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
3205 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
3206
3207             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3208      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
3209 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3210 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
3211 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3212 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3213 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
3214 c   Instead of res sim other local measure of b/b str reliability possible
3215             if (sigma_dih(k,i).ne.0)
3216      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3217 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3218           enddo
3219           lim_dih=nct-nnt-2
3220         endif
3221
3222         if (waga_theta.gt.0.0d0) then
3223 c         open (ientin,file=tpl_k_sigma_theta,status='old')
3224 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3225 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3226 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3227 c    &                              sigma_theta(k,i+nnt-1)
3228 c         enddo
3229 c1403   continue
3230 c         close (ientin)
3231
3232           do i = nnt+2,nct ! right? without parallel.
3233 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
3234 c         do i=ithet_start,ithet_end ! with FG parallel.
3235              if (idomain(k,i).eq.0) then
3236               sigma_theta(k,i)=0.0
3237               cycle
3238              endif
3239              thetatpl(k,i)=thetaref(i)
3240 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3241 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
3242 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
3243 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
3244 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
3245              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3246      &                        rescore(k,i-2))/3.0
3247 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3248              if (sigma_theta(k,i).ne.0)
3249      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3250
3251 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3252 c                             rescore(k,i-2) !  right expression ?
3253 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3254           enddo
3255         endif
3256
3257         if (waga_d.gt.0.0d0) then
3258 c       open (ientin,file=tpl_k_sigma_d,status='old')
3259 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3260 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3261 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3262 c    &                          sigma_d(k,i+nnt-1)
3263 c         enddo
3264 c1404   continue
3265
3266           do i = nnt,nct ! right? without parallel.
3267 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
3268 c         do i=loc_start,loc_end ! with FG parallel.
3269                if (itype(i).eq.10) cycle
3270                if (idomain(k,i).eq.0 ) then
3271                   sigma_d(k,i)=0.0
3272                   cycle
3273                endif
3274                xxtpl(k,i)=xxref(i)
3275                yytpl(k,i)=yyref(i)
3276                zztpl(k,i)=zzref(i)
3277 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3278 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3279 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3280 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
3281                sigma_d(k,i)=rescore3(k,i) !  right expression ?
3282                if (sigma_d(k,i).ne.0)
3283      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3284
3285 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
3286 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3287 c              read (ientin,*) sigma_d(k,i) ! 1st variant
3288           enddo
3289         endif
3290       enddo
3291 c
3292 c remove distance restraints not used in any model from the list
3293 c shift data in all arrays
3294 c
3295       if (waga_dist.ne.0.0d0) then
3296         ii=0
3297         liiflag=.true.
3298         do i=nnt,nct-2
3299          do j=i+2,nct
3300           ii=ii+1
3301           if (ii_in_use(ii).eq.0.and.liiflag) then
3302             liiflag=.false.
3303             iistart=ii
3304           endif
3305           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3306      &                   .not.liiflag.and.ii.eq.lim_odl) then
3307              if (ii.eq.lim_odl) then
3308               iishift=ii-iistart+1
3309              else
3310               iishift=ii-iistart
3311              endif
3312              liiflag=.true.
3313              do ki=iistart,lim_odl-iishift
3314               ires_homo(ki)=ires_homo(ki+iishift)
3315               jres_homo(ki)=jres_homo(ki+iishift)
3316               ii_in_use(ki)=ii_in_use(ki+iishift)
3317               do k=1,constr_homology
3318                odl(k,ki)=odl(k,ki+iishift)
3319                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3320                l_homo(k,ki)=l_homo(k,ki+iishift)
3321               enddo
3322              enddo
3323              ii=ii-iishift
3324              lim_odl=lim_odl-iishift
3325           endif
3326          enddo
3327         enddo
3328       endif
3329
3330       endif ! .not. klapaucjusz
3331
3332       if (constr_homology.gt.0) call homology_partition
3333       if (constr_homology.gt.0) call init_int_table
3334 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3335 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3336 c
3337 c Print restraints
3338 c
3339       if (.not.out_template_restr) return
3340 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3341       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3342        write (iout,*) "Distance restraints from templates"
3343        do ii=1,lim_odl
3344        write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3345      &  ii,ires_homo(ii),jres_homo(ii),
3346      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3347      &  ki=1,constr_homology)
3348        enddo
3349        write (iout,*) "Dihedral angle restraints from templates"
3350        do i=nnt+3,nct
3351         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3352      &      (rad2deg*dih(ki,i),
3353      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3354        enddo
3355        write (iout,*) "Virtual-bond angle restraints from templates"
3356        do i=nnt+2,nct
3357         write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3358      &      (rad2deg*thetatpl(ki,i),
3359      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3360        enddo
3361        write (iout,*) "SC restraints from templates"
3362        do i=nnt,nct
3363         write(iout,'(i5,100(4f8.2,4x))') i,
3364      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3365      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3366        enddo
3367       endif
3368 c -----------------------------------------------------------------
3369       return
3370       end
3371 c----------------------------------------------------------------------
3372 #ifdef WINIFL
3373       subroutine flush(iu)
3374       return
3375       end
3376 #endif
3377 #ifdef AIX
3378       subroutine flush(iu)
3379       call flush_(iu)
3380       return
3381       end
3382 #endif
3383 c------------------------------------------------------------------------------
3384       subroutine copy_to_tmp(source)
3385       implicit none
3386       include "DIMENSIONS"
3387       include "COMMON.IOUNITS"
3388       character*(*) source
3389       character* 256 tmpfile
3390       integer ilen
3391       external ilen
3392       logical ex
3393       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3394       inquire(file=tmpfile,exist=ex)
3395       if (ex) then
3396         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3397      &   " to temporary directory..."
3398         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3399         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3400       endif
3401       return
3402       end
3403 c------------------------------------------------------------------------------
3404       subroutine move_from_tmp(source)
3405       implicit none
3406       include "DIMENSIONS"
3407       include "COMMON.IOUNITS"
3408       character*(*) source
3409       integer ilen
3410       external ilen
3411       write (*,*) "Moving ",source(:ilen(source)),
3412      & " from temporary directory to working directory"
3413       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3414       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3415       return
3416       end
3417 c------------------------------------------------------------------------------
3418       subroutine random_init(seed)
3419 C
3420 C Initialize random number generator
3421 C
3422       implicit none
3423       include 'DIMENSIONS'
3424 #ifdef MPI
3425       include 'mpif.h'
3426       logical OKRandom, prng_restart
3427       real*8  r1
3428       integer iseed_array(4)
3429       integer error_msg,ierr
3430 #endif
3431       include 'COMMON.IOUNITS'
3432       include 'COMMON.TIME1'
3433       include 'COMMON.THREAD'
3434       include 'COMMON.SBRIDGE'
3435       include 'COMMON.CONTROL'
3436       include 'COMMON.MCM'
3437       include 'COMMON.MAP'
3438       include 'COMMON.HEADER'
3439       include 'COMMON.CSA'
3440       include 'COMMON.CHAIN'
3441       include 'COMMON.MUCA'
3442       include 'COMMON.MD'
3443       include 'COMMON.FFIELD'
3444       include 'COMMON.SETUP'
3445       integer i,iseed
3446       double precision seed,ran_number
3447       iseed=-dint(dabs(seed))
3448       if (iseed.eq.0) then
3449         write (iout,'(/80(1h*)/20x,a/80(1h*))')
3450      &    'Random seed undefined. The program will stop.'
3451         write (*,'(/80(1h*)/20x,a/80(1h*))')
3452      &    'Random seed undefined. The program will stop.'
3453 #ifdef MPI
3454         call mpi_finalize(mpi_comm_world,ierr)
3455 #endif
3456         stop 'Bad random seed.'
3457       endif
3458 #ifdef MPI
3459       if (fg_rank.eq.0) then
3460       seed=seed*(me+1)+1
3461 #ifdef AMD64
3462       if(me.eq.king)
3463      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3464       OKRandom = prng_restart(me,iseed)
3465 #else
3466       do i=1,4
3467        tmp=65536.0d0**(4-i)
3468        iseed_array(i) = dint(seed/tmp)
3469        seed=seed-iseed_array(i)*tmp
3470       enddo
3471       if(me.eq.king)
3472      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3473      &                 (iseed_array(i),i=1,4)
3474       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3475      &                 (iseed_array(i),i=1,4)
3476       OKRandom = prng_restart(me,iseed_array)
3477 #endif
3478       if (OKRandom) then
3479 c        r1 = prng_next(me)
3480         r1=ran_number(0.0D0,1.0D0)
3481         if(me.eq.king)
3482      &   write (iout,*) 'ran_num',r1
3483         if (r1.lt.0.0d0) OKRandom=.false.
3484       endif
3485       if (.not.OKRandom) then
3486         write (iout,*) 'PRNG IS NOT WORKING!!!'
3487         print *,'PRNG IS NOT WORKING!!!'
3488         if (me.eq.0) then
3489          call flush(iout)
3490          call mpi_abort(mpi_comm_world,error_msg,ierr)
3491          stop
3492         else
3493          write (iout,*) 'too many processors for parallel prng'
3494          write (*,*) 'too many processors for parallel prng'
3495          call flush(iout)
3496          stop
3497         endif
3498       endif
3499       endif
3500 #else
3501       call vrndst(iseed)
3502       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3503 #endif
3504       return
3505       end
3506 c----------------------------------------------------------------------
3507       subroutine read_klapaucjusz
3508       implicit none
3509       include 'DIMENSIONS'
3510 #ifdef MPI
3511       include 'mpif.h'
3512 #endif
3513       include 'COMMON.SETUP'
3514       include 'COMMON.CONTROL'
3515       include 'COMMON.HOMOLOGY'
3516       include 'COMMON.CHAIN'
3517       include 'COMMON.IOUNITS'
3518       include 'COMMON.MD'
3519       include 'COMMON.GEO'
3520       include 'COMMON.INTERACT'
3521       include 'COMMON.NAMES'
3522       character*256 fragfile
3523       integer ninclust(maxclust),inclust(max_template,maxclust),
3524      &  nresclust(maxclust),iresclust(maxres,maxclust),nclust
3525
3526       character*2 kic2
3527       character*24 model_ki_dist, model_ki_angle
3528       character*500 controlcard
3529       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp,
3530      & ik,ll,ii,kk,iistart,iishift,lim_xx
3531       double precision distal
3532       logical lprn /.true./
3533       integer ilen
3534       external ilen
3535       logical liiflag
3536 c
3537 c
3538       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3539       double precision, dimension (max_template,maxres) :: rescore
3540       double precision, dimension (max_template,maxres) :: rescore2
3541       character*24 pdbfile,tpl_k_rescore
3542
3543 c
3544 c For new homol impl
3545 c
3546       include 'COMMON.VAR'
3547 c
3548       call getenv("FRAGFILE",fragfile)
3549       open(ientin,file=fragfile,status="old",err=10)
3550       read(ientin,*) constr_homology,nclust
3551       l_homo = .false.
3552       sigma_theta=0.0
3553       sigma_d=0.0
3554       sigma_dih=0.0
3555 c Read pdb files
3556       do k=1,constr_homology
3557         read(ientin,'(a)') pdbfile
3558         write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3559      &  pdbfile(:ilen(pdbfile))
3560         open(ipdbin,file=pdbfile,status='old',err=33)
3561         goto 34
3562   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
3563      &  pdbfile(:ilen(pdbfile))
3564         stop
3565   34    continue
3566         unres_pdb=.false.
3567         call readpdb_template(k)
3568         do i=1,nres
3569           rescore(k,i)=0.2d0
3570           rescore2(k,i)=1.0d0
3571         enddo
3572       enddo
3573 c Read clusters
3574       do i=1,nclust
3575         read(ientin,*) ninclust(i),nresclust(i)
3576         read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3577         read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3578       enddo
3579 c
3580 c Loop over clusters
3581 c
3582       do l=1,nclust
3583         do ll = 1,ninclust(l)
3584
3585         k = inclust(ll,l)
3586         do i=1,nres
3587           idomain(k,i)=0
3588         enddo
3589         do i=1,nresclust(l)
3590           if (nnt.gt.1)  then
3591             idomain(k,iresclust(i,l)+1) = 1
3592           else
3593             idomain(k,iresclust(i,l)) = 1
3594           endif
3595         enddo
3596 c
3597 c     Distance restraints
3598 c
3599 c          ... --> odl(k,ii)
3600 C Copy the coordinates from reference coordinates (?)
3601         do i=1,2*nres
3602           do j=1,3
3603             c(j,i)=chomo(j,i,k)
3604 c           write (iout,*) "c(",j,i,") =",c(j,i)
3605           enddo
3606         enddo
3607         call int_from_cart(.true.,.false.)
3608         call sc_loc_geom(.true.)
3609         do i=1,nres
3610           thetaref(i)=theta(i)
3611           phiref(i)=phi(i)
3612         enddo
3613         if (waga_dist.ne.0.0d0) then
3614           ii=0
3615           do i = nnt,nct-2
3616             do j=i+2,nct
3617
3618               x12=c(1,i)-c(1,j)
3619               y12=c(2,i)-c(2,j)
3620               z12=c(3,i)-c(3,j)
3621               distal=dsqrt(x12*x12+y12*y12+z12*z12)
3622 c              write (iout,*) k,i,j,distal,dist2_cut
3623
3624             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3625      &            .and. distal.le.dist2_cut ) then
3626
3627               ii=ii+1
3628               ii_in_use(ii)=1
3629               l_homo(k,ii)=.true.
3630
3631 c             write (iout,*) "k",k
3632 c             write (iout,*) "i",i," j",j," constr_homology",
3633 c    &                       constr_homology
3634               ires_homo(ii)=i
3635               jres_homo(ii)=j
3636               odl(k,ii)=distal
3637               if (read2sigma) then
3638                 sigma_odl(k,ii)=0
3639                 do ik=i,j
3640                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3641                 enddo
3642                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3643                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3644      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3645               else
3646                 if (odl(k,ii).le.dist_cut) then
3647                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3648                 else
3649 #ifdef OLDSIGMA
3650                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3651      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3652 #else
3653                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3654      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3655 #endif
3656                 endif
3657               endif
3658               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3659             else
3660               ii=ii+1
3661 c              l_homo(k,ii)=.false.
3662             endif
3663             enddo
3664           enddo
3665         lim_odl=ii
3666         endif
3667 c
3668 c     Theta, dihedral and SC retraints
3669 c
3670         if (waga_angle.gt.0.0d0) then
3671           do i = nnt+3,nct
3672             if (idomain(k,i).eq.0) then
3673 c               sigma_dih(k,i)=0.0
3674                cycle
3675             endif
3676             dih(k,i)=phiref(i)
3677             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3678      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
3679 c            write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3680 c     &       " sigma_dihed",sigma_dih(k,i)
3681             if (sigma_dih(k,i).ne.0)
3682      &       sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3683           enddo
3684           lim_dih=nct-nnt-2
3685         endif
3686
3687         if (waga_theta.gt.0.0d0) then
3688           do i = nnt+2,nct
3689              if (idomain(k,i).eq.0) then
3690 c              sigma_theta(k,i)=0.0
3691               cycle
3692              endif
3693              thetatpl(k,i)=thetaref(i)
3694              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3695      &                        rescore(k,i-2))/3.0
3696              if (sigma_theta(k,i).ne.0)
3697      &       sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3698           enddo
3699         endif
3700
3701         if (waga_d.gt.0.0d0) then
3702           do i = nnt,nct
3703                if (itype(i).eq.10) cycle
3704                if (idomain(k,i).eq.0 ) then
3705 c                  sigma_d(k,i)=0.0
3706                   cycle
3707                endif
3708                xxtpl(k,i)=xxref(i)
3709                yytpl(k,i)=yyref(i)
3710                zztpl(k,i)=zzref(i)
3711                sigma_d(k,i)=rescore(k,i)
3712                if (sigma_d(k,i).ne.0)
3713      &          sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3714                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3715           enddo
3716         endif
3717       enddo ! l
3718       enddo ! ll
3719 c
3720 c remove distance restraints not used in any model from the list
3721 c shift data in all arrays
3722 c
3723       if (waga_dist.ne.0.0d0) then
3724         ii=0
3725         liiflag=.true.
3726         do i=nnt,nct-2
3727          do j=i+2,nct
3728           ii=ii+1
3729           if (ii_in_use(ii).eq.0.and.liiflag) then
3730             liiflag=.false.
3731             iistart=ii
3732           endif
3733           if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3734      &                   .not.liiflag.and.ii.eq.lim_odl) then
3735              if (ii.eq.lim_odl) then
3736               iishift=ii-iistart+1
3737              else
3738               iishift=ii-iistart
3739              endif
3740              liiflag=.true.
3741              do ki=iistart,lim_odl-iishift
3742               ires_homo(ki)=ires_homo(ki+iishift)
3743               jres_homo(ki)=jres_homo(ki+iishift)
3744               ii_in_use(ki)=ii_in_use(ki+iishift)
3745               do k=1,constr_homology
3746                odl(k,ki)=odl(k,ki+iishift)
3747                sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3748                l_homo(k,ki)=l_homo(k,ki+iishift)
3749               enddo
3750              enddo
3751              ii=ii-iishift
3752              lim_odl=lim_odl-iishift
3753           endif
3754          enddo
3755         enddo
3756       endif
3757
3758       return
3759    10 stop "Error in fragment file"
3760       end