source/unres/src_MD-M/readrtns_CSA.F

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