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