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)
  45      &write (iout,'(a,i5,a)') "The following",nhpb-nss,
  46      & " distance constraints have been imposed"
  47       do i=nss+1,nhpb
  48         write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
  49       enddo
  50 #ifdef MPI
  51       endif
  52 #endif
  53 c      print *,"Processor",myrank," leaves READRTNS"
  54       return
  55       end
  56 C-------------------------------------------------------------------------------
  57       subroutine read_control
  58 C
  59 C Read contorl data
  60 C
  61       implicit real*8 (a-h,o-z)
  62       include 'DIMENSIONS'
  63 #ifdef MP
  64       include 'mpif.h'
  65       logical OKRandom, prng_restart
  66       real*8  r1
  67 #endif
  68       include 'COMMON.IOUNITS'
  69       include 'COMMON.TIME1'
  70       include 'COMMON.THREAD'
  71       include 'COMMON.SBRIDGE'
  72       include 'COMMON.CONTROL'
  73       include 'COMMON.MCM'
  74       include 'COMMON.MAP'
  75       include 'COMMON.HEADER'
  76       include 'COMMON.CSA'
  77       include 'COMMON.CHAIN'
  78       include 'COMMON.MUCA'
  79       include 'COMMON.MD'
  80       include 'COMMON.FFIELD'
  81       include 'COMMON.INTERACT'
  82       include 'COMMON.SETUP'
  83       include 'COMMON.SPLITELE'
  84       include 'COMMON.SHIELD'
  85       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
  86       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
  87       character*80 ucase
  88       character*320 controlcard
  89
  90       nglob_csa=0
  91       eglob_csa=1d99
  92       nmin_csa=0
  93       read (INP,'(a)') titel
  94       call card_concat(controlcard)
  95 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
  96 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
  97       call reada(controlcard,'SEED',seed,0.0D0)
  98       call random_init(seed)
  99 C Set up the time limit (caution! The time must be input in minutes!)
 100       read_cart=index(controlcard,'READ_CART').gt.0
 101       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
 102 C this variable with_theta_constr is the variable which allow to read and execute the
 103 C constrains on theta angles WITH_THETA_CONSTR is the keyword
 104       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
 105       write (iout,*) "with_theta_constr ",with_theta_constr
 106       call readi(controlcard,'SYM',symetr,1)
 107       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
 108       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
 109       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
 110       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
 111       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
 112       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
 113       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 114       call reada(controlcard,'DRMS',drms,0.1D0)
 115       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 116        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 117        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 118        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 119        write (iout,'(a,f10.1)')'DRMS    = ',drms
 120        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 121        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 122       endif
 123       call readi(controlcard,'NZ_START',nz_start,0)
 124       call readi(controlcard,'NZ_END',nz_end,0)
 125       call readi(controlcard,'IZ_SC',iz_sc,0)
 126       timlim=60.0D0*timlim
 127       safety = 60.0d0*safety
 128       timem=timlim
 129       modecalc=0
 130       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 131       minim=(index(controlcard,'MINIMIZE').gt.0)
 132       dccart=(index(controlcard,'CART').gt.0)
 133       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 134       overlapsc=.not.overlapsc
 135       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
 136       searchsc=.not.searchsc
 137       sideadd=(index(controlcard,'SIDEADD').gt.0)
 138       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 139       outpdb=(index(controlcard,'PDBOUT').gt.0)
 140       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 141       pdbref=(index(controlcard,'PDBREF').gt.0)
 142       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 143       indpdb=index(controlcard,'PDBSTART')
 144       extconf=(index(controlcard,'EXTCONF').gt.0)
 145       AFMlog=(index(controlcard,'AFM'))
 146       selfguide=(index(controlcard,'SELFGUIDE'))
 147       print *,'AFMlog',AFMlog,selfguide,"KUPA"
 148       call readi(controlcard,'IPRINT',iprint,0)
 149 C SHIELD keyword sets if the shielding effect of side-chains is used
 150 C 0 denots no shielding is used all peptide are equally despite the
 151 C solvent accesible area
 152 C 1 the newly introduced function
 153 C 2 reseved for further possible developement
 154       call readi(controlcard,'SHIELD',shield_mode,0)
 155 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 156         write(iout,*) "shield_mode",shield_mode
 157 C      endif
 158       call readi(controlcard,'MAXGEN',maxgen,10000)
 159       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 160       call readi(controlcard,"KDIAG",kdiag,0)
 161       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
 162       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 163      & write (iout,*) "RESCALE_MODE",rescale_mode
 164       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 165       if (index(controlcard,'REGULAR').gt.0.0D0) then
 166         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 167         modecalc=1
 168         refstr=.true.
 169       endif
 170       if (index(controlcard,'CHECKGRAD').gt.0) then
 171         modecalc=5
 172         if (index(controlcard,'CART').gt.0) then
 173           icheckgrad=1
 174         elseif (index(controlcard,'CARINT').gt.0) then
 175           icheckgrad=2
 176         else
 177           icheckgrad=3
 178         endif
 179       elseif (index(controlcard,'THREAD').gt.0) then
 180         modecalc=2
 181         call readi(controlcard,'THREAD',nthread,0)
 182         if (nthread.gt.0) then
 183           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 184         else
 185           if (fg_rank.eq.0)
 186      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 187           stop 'Error termination in Read_Control.'
 188         endif
 189       else if (index(controlcard,'MCMA').gt.0) then
 190         modecalc=3
 191       else if (index(controlcard,'MCEE').gt.0) then
 192         modecalc=6
 193       else if (index(controlcard,'MULTCONF').gt.0) then
 194         modecalc=4
 195       else if (index(controlcard,'MAP').gt.0) then
 196         modecalc=7
 197         call readi(controlcard,'MAP',nmap,0)
 198       else if (index(controlcard,'CSA').gt.0) then
 199         modecalc=8
 200 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 201 crc
 202 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 203 crc
 204 crc        modecalc=9
 205 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 206 cfmc        modecalc=10
 207       else if (index(controlcard,'SOFTREG').gt.0) then
 208         modecalc=11
 209       else if (index(controlcard,'CHECK_BOND').gt.0) then
 210         modecalc=-1
 211       else if (index(controlcard,'TEST').gt.0) then
 212         modecalc=-2
 213       else if (index(controlcard,'MD').gt.0) then
 214         modecalc=12
 215       else if (index(controlcard,'RE ').gt.0) then
 216         modecalc=14
 217       endif
 218
 219       lmuca=index(controlcard,'MUCA').gt.0
 220       call readi(controlcard,'MUCADYN',mucadyn,0)
 221       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 222       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 223      & then
 224        write (iout,*) 'MUCADYN=',mucadyn
 225        write (iout,*) 'MUCASMOOTH=',muca_smooth
 226       endif
 227
 228       iscode=index(controlcard,'ONE_LETTER')
 229       indphi=index(controlcard,'PHI')
 230       indback=index(controlcard,'BACK')
 231       iranconf=index(controlcard,'RAND_CONF')
 232       i2ndstr=index(controlcard,'USE_SEC_PRED')
 233       gradout=index(controlcard,'GRADOUT').gt.0
 234       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 235 C  DISTCHAINMAX become obsolete for periodic boundry condition
 236       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
 237 C Reading the dimensions of box in x,y,z coordinates
 238       call reada(controlcard,'BOXX',boxxsize,100.0d0)
 239       call reada(controlcard,'BOXY',boxysize,100.0d0)
 240       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
 241 c Cutoff range for interactions
 242       call reada(controlcard,"R_CUT",r_cut,15.0d0)
 243       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 244       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
 245       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
 246       if (lipthick.gt.0.0d0) then
 247        bordliptop=(boxzsize+lipthick)/2.0
 248        bordlipbot=bordliptop-lipthick
 249 C      endif
 250       if ((bordliptop.gt.boxzsize).or.(borlipbot.lt.0.0))
 251      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
 252       buflipbot=bordlipbot+lipbufthick
 253       bufliptop=bordliptop-lipbufthick
 254       if ((lipbufthick*2.0d0).gt.lipthick)
 255      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
 256       endif
 257       write(iout,*) "bordliptop=",bordliptop
 258       write(iout,*) "bordlipbot=",bordlipbot
 259       write(iout,*) "bufliptop=",bufliptop
 260       write(iout,*) "buflipbot=",buflipbot
 261       write (iout,*) "SHIELD MODE",shield_mode
 262       if (shield_mode.gt.0) then
 263       pi=3.141592d0
 264 C VSolvSphere the volume of solving sphere
 265 C      print *,pi,"pi"
 266 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
 267 C there will be no distinction between proline peptide group and normal peptide
 268 C group in case of shielding parameters
 269       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
 270       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
 271       write (iout,*) VSolvSphere,VSolvSphere_div
 272 C long axis of side chain
 273       do i=1,ntyp
 274       long_r_sidechain(i)=vbldsc0(1,i)
 275       short_r_sidechain(i)=sigma0(i)
 276       enddo
 277       buff_shield=1.0d0
 278       endif
 279       if (me.eq.king .or. .not.out1file )
 280      &  write (iout,*) "DISTCHAINMAX",distchainmax
 281
 282       if(me.eq.king.or..not.out1file)
 283      & write (iout,'(2a)') diagmeth(kdiag),
 284      &  ' routine used to diagonalize matrices.'
 285       return
 286       end
 287 c--------------------------------------------------------------------------
 288       subroutine read_REMDpar
 289 C
 290 C Read REMD settings
 291 C
 292       implicit real*8 (a-h,o-z)
 293       include 'DIMENSIONS'
 294       include 'COMMON.IOUNITS'
 295       include 'COMMON.TIME1'
 296       include 'COMMON.MD'
 297 #ifndef LANG0
 298       include 'COMMON.LANGEVIN'
 299 #else
 300       include 'COMMON.LANGEVIN.lang0'
 301 #endif
 302       include 'COMMON.INTERACT'
 303       include 'COMMON.NAMES'
 304       include 'COMMON.GEO'
 305       include 'COMMON.REMD'
 306       include 'COMMON.CONTROL'
 307       include 'COMMON.SETUP'
 308       character*80 ucase
 309       character*320 controlcard
 310       character*3200 controlcard1
 311       integer iremd_m_total
 312
 313       if(me.eq.king.or..not.out1file)
 314      & write (iout,*) "REMD setup"
 315
 316       call card_concat(controlcard)
 317       call readi(controlcard,"NREP",nrep,3)
 318       call readi(controlcard,"NSTEX",nstex,1000)
 319       call reada(controlcard,"RETMIN",retmin,10.0d0)
 320       call reada(controlcard,"RETMAX",retmax,1000.0d0)
 321       mremdsync=(index(controlcard,'SYNC').gt.0)
 322       call readi(controlcard,"NSYN",i_sync_step,100)
 323       restart1file=(index(controlcard,'REST1FILE').gt.0)
 324       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
 325       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
 326       if(max_cache_traj_use.gt.max_cache_traj)
 327      &           max_cache_traj_use=max_cache_traj
 328       if(me.eq.king.or..not.out1file) then
 329 cd       if (traj1file) then
 330 crc caching is in testing - NTWX is not ignored
 331 cd        write (iout,*) "NTWX value is ignored"
 332 cd        write (iout,*) "  trajectory is stored to one file by master"
 333 cd        write (iout,*) "  before exchange at NSTEX intervals"
 334 cd       endif
 335        write (iout,*) "NREP= ",nrep
 336        write (iout,*) "NSTEX= ",nstex
 337        write (iout,*) "SYNC= ",mremdsync
 338        write (iout,*) "NSYN= ",i_sync_step
 339        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
 340       endif
 341       remd_tlist=.false.
 342       if (index(controlcard,'TLIST').gt.0) then
 343          remd_tlist=.true.
 344          call card_concat(controlcard1)
 345          read(controlcard1,*) (remd_t(i),i=1,nrep)
 346          if(me.eq.king.or..not.out1file)
 347      &    write (iout,*)'tlist',(remd_t(i),i=1,nrep)
 348       endif
 349       remd_mlist=.false.
 350       if (index(controlcard,'MLIST').gt.0) then
 351          remd_mlist=.true.
 352          call card_concat(controlcard1)
 353          read(controlcard1,*) (remd_m(i),i=1,nrep)
 354          if(me.eq.king.or..not.out1file) then
 355           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
 356           iremd_m_total=0
 357           do i=1,nrep
 358            iremd_m_total=iremd_m_total+remd_m(i)
 359           enddo
 360            write (iout,*) 'Total number of replicas ',iremd_m_total
 361           endif
 362          endif
 363       if(me.eq.king.or..not.out1file)
 364      &   write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
 365       return
 366       end
 367 c--------------------------------------------------------------------------
 368       subroutine read_MDpar
 369 C
 370 C Read MD settings
 371 C
 372       implicit real*8 (a-h,o-z)
 373       include 'DIMENSIONS'
 374       include 'COMMON.IOUNITS'
 375       include 'COMMON.TIME1'
 376       include 'COMMON.MD'
 377 #ifndef LANG0
 378       include 'COMMON.LANGEVIN'
 379 #else
 380       include 'COMMON.LANGEVIN.lang0'
 381 #endif
 382       include 'COMMON.INTERACT'
 383       include 'COMMON.NAMES'
 384       include 'COMMON.GEO'
 385       include 'COMMON.SETUP'
 386       include 'COMMON.CONTROL'
 387       include 'COMMON.SPLITELE'
 388       character*80 ucase
 389       character*320 controlcard
 390
 391       call card_concat(controlcard)
 392       call readi(controlcard,"NSTEP",n_timestep,1000000)
 393       call readi(controlcard,"NTWE",ntwe,100)
 394       call readi(controlcard,"NTWX",ntwx,1000)
 395       call reada(controlcard,"DT",d_time,1.0d-1)
 396       call reada(controlcard,"DVMAX",dvmax,2.0d1)
 397       call reada(controlcard,"DAMAX",damax,1.0d1)
 398       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
 399       call readi(controlcard,"LANG",lang,0)
 400       RESPA = index(controlcard,"RESPA") .gt. 0
 401       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
 402       ntime_split0=ntime_split
 403       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
 404       ntime_split0=ntime_split
 405 c      call reada(controlcard,"R_CUT",r_cut,2.0d0)
 406 c      call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 407       rest = index(controlcard,"REST").gt.0
 408       tbf = index(controlcard,"TBF").gt.0
 409       usampl = index(controlcard,"USAMPL").gt.0
 410
 411       mdpdb = index(controlcard,"MDPDB").gt.0
 412       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 413       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
 414       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
 415       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
 416       if (count_reset_moment.eq.0) count_reset_moment=1000000000
 417       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
 418       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
 419       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
 420       if (count_reset_vel.eq.0) count_reset_vel=1000000000
 421       large = index(controlcard,"LARGE").gt.0
 422       print_compon = index(controlcard,"PRINT_COMPON").gt.0
 423       rattle = index(controlcard,"RATTLE").gt.0
 424 c  if performing umbrella sampling, fragments constrained are read from the fragment file
 425       nset=0
 426       if(usampl) then
 427         call read_fragments
 428       endif
 429
 430       if(me.eq.king.or..not.out1file) then
 431        write (iout,*)
 432        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
 433        write (iout,*)
 434        write (iout,'(a)') "The units are:"
 435        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
 436        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
 437      &  " acceleration: angstrom/(48.9 fs)**2"
 438        write (iout,'(a)') "energy: kcal/mol, temperature: K"
 439        write (iout,*)
 440        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
 441        write (iout,'(a60,f10.5,a)')
 442      &  "Initial time step of numerical integration:",d_time,
 443      &  " natural units"
 444        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
 445        if (RESPA) then
 446         write (iout,'(2a,i4,a)')
 447      &    "A-MTS algorithm used; initial time step for fast-varying",
 448      &    " short-range forces split into",ntime_split," steps."
 449         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
 450      &   r_cut," lambda",rlamb
 451        endif
 452        write (iout,'(2a,f10.5)')
 453      &  "Maximum acceleration threshold to reduce the time step",
 454      &  "/increase split number:",damax
 455        write (iout,'(2a,f10.5)')
 456      &  "Maximum predicted energy drift to reduce the timestep",
 457      &  "/increase split number:",edriftmax
 458        write (iout,'(a60,f10.5)')
 459      & "Maximum velocity threshold to reduce velocities:",dvmax
 460        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
 461        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
 462        if (rattle) write (iout,'(a60)')
 463      &  "Rattle algorithm used to constrain the virtual bonds"
 464       endif
 465       reset_fricmat=1000
 466       if (lang.gt.0) then
 467         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
 468         call reada(controlcard,"RWAT",rwat,1.4d0)
 469         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
 470         surfarea=index(controlcard,"SURFAREA").gt.0
 471         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
 472         if(me.eq.king.or..not.out1file)then
 473          write (iout,'(/a,$)') "Langevin dynamics calculation"
 474          if (lang.eq.1) then
 475           write (iout,'(a/)')
 476      &      " with direct integration of Langevin equations"
 477          else if (lang.eq.2) then
 478           write (iout,'(a/)') " with TINKER stochasic MD integrator"
 479          else if (lang.eq.3) then
 480           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
 481          else if (lang.eq.4) then
 482           write (iout,'(a/)') " in overdamped mode"
 483          else
 484           write (iout,'(//a,i5)')
 485      &      "=========== ERROR: Unknown Langevin dynamics mode:",lang
 486           stop
 487          endif
 488          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 489          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
 490          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
 491          write (iout,'(a60,f10.5)')
 492      &   "Scaling factor of the friction forces:",scal_fric
 493          if (surfarea) write (iout,'(2a,i10,a)')
 494      &     "Friction coefficients will be scaled by solvent-accessible",
 495      &     " surface area every",reset_fricmat," steps."
 496         endif
 497 c Calculate friction coefficients and bounds of stochastic forces
 498         eta=6*pi*cPoise*etawat
 499         if(me.eq.king.or..not.out1file)
 500      &   write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
 501      &   ,eta
 502         gamp=scal_fric*(pstok+rwat)*eta
 503         stdfp=dsqrt(2*Rb*t_bath/d_time)
 504         do i=1,ntyp
 505           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
 506           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 507         enddo
 508         if(me.eq.king.or..not.out1file)then
 509          write (iout,'(/2a/)')
 510      &   "Radii of site types and friction coefficients and std's of",
 511      &   " stochastic forces of fully exposed sites"
 512          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
 513          do i=1,ntyp
 514           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
 515      &     gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
 516          enddo
 517         endif
 518       else if (tbf) then
 519         if(me.eq.king.or..not.out1file)then
 520          write (iout,'(a)') "Berendsen bath calculation"
 521          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 522          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
 523          if (reset_moment)
 524      &   write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
 525      &   count_reset_moment," steps"
 526          if (reset_vel)
 527      &    write (iout,'(a,i10,a)')
 528      &    "Velocities will be reset at random every",count_reset_vel,
 529      &   " steps"
 530         endif
 531       else
 532         if(me.eq.king.or..not.out1file)
 533      &   write (iout,'(a31)') "Microcanonical mode calculation"
 534       endif
 535       if(me.eq.king.or..not.out1file)then
 536        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
 537        if (usampl) then
 538           write(iout,*) "MD running with constraints."
 539           write(iout,*) "Equilibration time ", eq_time, " mtus."
 540           write(iout,*) "Constraining ", nfrag," fragments."
 541           write(iout,*) "Length of each fragment, weight and q0:"
 542           do iset=1,nset
 543            write (iout,*) "Set of restraints #",iset
 544            do i=1,nfrag
 545               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
 546      &           ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
 547            enddo
 548            write(iout,*) "constraints between ", npair, "fragments."
 549            write(iout,*) "constraint pairs, weights and q0:"
 550            do i=1,npair
 551             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
 552      &             ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
 553            enddo
 554            write(iout,*) "angle constraints within ", nfrag_back,
 555      &      "backbone fragments."
 556            write(iout,*) "fragment, weights:"
 557            do i=1,nfrag_back
 558             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
 559      &         ifrag_back(2,i,iset),wfrag_back(1,i,iset),
 560      &         wfrag_back(2,i,iset),wfrag_back(3,i,iset)
 561            enddo
 562           enddo
 563         iset=mod(kolor,nset)+1
 564        endif
 565       endif
 566       if(me.eq.king.or..not.out1file)
 567      & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
 568       return
 569       end
 570 c------------------------------------------------------------------------------
 571       subroutine molread
 572 C
 573 C Read molecular data.
 574 C
 575       implicit real*8 (a-h,o-z)
 576       include 'DIMENSIONS'
 577 #ifdef MPI
 578       include 'mpif.h'
 579       integer error_msg
 580 #endif
 581       include 'COMMON.IOUNITS'
 582       include 'COMMON.GEO'
 583       include 'COMMON.VAR'
 584       include 'COMMON.INTERACT'
 585       include 'COMMON.LOCAL'
 586       include 'COMMON.NAMES'
 587       include 'COMMON.CHAIN'
 588       include 'COMMON.FFIELD'
 589       include 'COMMON.SBRIDGE'
 590       include 'COMMON.HEADER'
 591       include 'COMMON.CONTROL'
 592       include 'COMMON.DBASE'
 593       include 'COMMON.THREAD'
 594       include 'COMMON.CONTACTS'
 595       include 'COMMON.TORCNSTR'
 596       include 'COMMON.TIME1'
 597       include 'COMMON.BOUNDS'
 598       include 'COMMON.MD'
 599       include 'COMMON.SETUP'
 600       character*4 sequence(maxres)
 601       integer rescode
 602       double precision x(maxvar)
 603       character*256 pdbfile
 604       character*400 weightcard
 605       character*80 weightcard_t,ucase
 606       dimension itype_pdb(maxres)
 607       common /pizda/ itype_pdb
 608       logical seq_comp,fail
 609       double precision energia(0:n_ene)
 610       integer ilen
 611       external ilen
 612 C
 613 C Body
 614 C
 615 C Read weights of the subsequent energy terms.
 616        call card_concat(weightcard)
 617        call reada(weightcard,'WLONG',wlong,1.0D0)
 618        call reada(weightcard,'WSC',wsc,wlong)
 619        call reada(weightcard,'WSCP',wscp,wlong)
 620        call reada(weightcard,'WELEC',welec,1.0D0)
 621        call reada(weightcard,'WVDWPP',wvdwpp,welec)
 622        call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 623        call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 624        call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 625        call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 626        call reada(weightcard,'WTURN3',wturn3,1.0D0)
 627        call reada(weightcard,'WTURN4',wturn4,1.0D0)
 628        call reada(weightcard,'WTURN6',wturn6,1.0D0)
 629        call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 630        call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 631        call reada(weightcard,'WBOND',wbond,1.0D0)
 632        call reada(weightcard,'WTOR',wtor,1.0D0)
 633        call reada(weightcard,'WTORD',wtor_d,1.0D0)
 634        call reada(weightcard,'WANG',wang,1.0D0)
 635        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 636        call reada(weightcard,'SCAL14',scal14,0.4D0)
 637        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 638        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 639        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 640        call reada(weightcard,'TEMP0',temp0,300.0d0)
 641        call reada(weightcard,'WLT',wliptran,0.0D0)
 642        if (index(weightcard,'SOFT').gt.0) ipot=6
 643 C 12/1/95 Added weight for the multi-body term WCORR
 644        call reada(weightcard,'WCORRH',wcorr,1.0D0)
 645        if (wcorr4.gt.0.0d0) wcorr=wcorr4
 646        weights(1)=wsc
 647        weights(2)=wscp
 648        weights(3)=welec
 649        weights(4)=wcorr
 650        weights(5)=wcorr5
 651        weights(6)=wcorr6
 652        weights(7)=wel_loc
 653        weights(8)=wturn3
 654        weights(9)=wturn4
 655        weights(10)=wturn6
 656        weights(11)=wang
 657        weights(12)=wscloc
 658        weights(13)=wtor
 659        weights(14)=wtor_d
 660        weights(15)=wstrain
 661        weights(16)=wvdwpp
 662        weights(17)=wbond
 663        weights(18)=scal14
 664        weights(21)=wsccor
 665       if(me.eq.king.or..not.out1file)
 666      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 667      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 668      &  wturn4,wturn6
 669    10 format (/'Energy-term weights (unscaled):'//
 670      & 'WSCC=   ',f10.6,' (SC-SC)'/
 671      & 'WSCP=   ',f10.6,' (SC-p)'/
 672      & 'WELEC=  ',f10.6,' (p-p electr)'/
 673      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 674      & 'WBOND=  ',f10.6,' (stretching)'/
 675      & 'WANG=   ',f10.6,' (bending)'/
 676      & 'WSCLOC= ',f10.6,' (SC local)'/
 677      & 'WTOR=   ',f10.6,' (torsional)'/
 678      & 'WTORD=  ',f10.6,' (double torsional)'/
 679      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 680      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 681      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 682      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 683      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 684      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
 685      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 686      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 687      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 688       if(me.eq.king.or..not.out1file)then
 689        if (wcorr4.gt.0.0d0) then
 690         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 691      &   'between contact pairs of peptide groups'
 692         write (iout,'(2(a,f5.3/))')
 693      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 694      &  'Range of quenching the correlation terms:',2*delt_corr
 695        else if (wcorr.gt.0.0d0) then
 696         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 697      &   'between contact pairs of peptide groups'
 698        endif
 699        write (iout,'(a,f8.3)')
 700      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 701        write (iout,'(a,f8.3)')
 702      &  'General scaling factor of SC-p interactions:',scalscp
 703       endif
 704       r0_corr=cutoff_corr-delt_corr
 705       do i=1,ntyp
 706         aad(i,1)=scalscp*aad(i,1)
 707         aad(i,2)=scalscp*aad(i,2)
 708         bad(i,1)=scalscp*bad(i,1)
 709         bad(i,2)=scalscp*bad(i,2)
 710       enddo
 711       call rescale_weights(t_bath)
 712       if(me.eq.king.or..not.out1file)
 713      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 714      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 715      &  wturn4,wturn6
 716    22 format (/'Energy-term weights (scaled):'//
 717      & 'WSCC=   ',f10.6,' (SC-SC)'/
 718      & 'WSCP=   ',f10.6,' (SC-p)'/
 719      & 'WELEC=  ',f10.6,' (p-p electr)'/
 720      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 721      & 'WBOND=  ',f10.6,' (stretching)'/
 722      & 'WANG=   ',f10.6,' (bending)'/
 723      & 'WSCLOC= ',f10.6,' (SC local)'/
 724      & 'WTOR=   ',f10.6,' (torsional)'/
 725      & 'WTORD=  ',f10.6,' (double torsional)'/
 726      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 727      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 728      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 729      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 730      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 731      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
 732      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 733      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 734      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 735       if(me.eq.king.or..not.out1file)
 736      & write (iout,*) "Reference temperature for weights calculation:",
 737      &  temp0
 738       call reada(weightcard,"D0CM",d0cm,3.78d0)
 739       call reada(weightcard,"AKCM",akcm,15.1d0)
 740       call reada(weightcard,"AKTH",akth,11.0d0)
 741       call reada(weightcard,"AKCT",akct,12.0d0)
 742       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 743       call reada(weightcard,"V2SS",v2ss,7.61d0)
 744       call reada(weightcard,"V3SS",v3ss,13.7d0)
 745       call reada(weightcard,"EBR",ebr,-5.50D0)
 746       call reada(weightcard,"ATRISS",atriss,0.301D0)
 747       call reada(weightcard,"BTRISS",btriss,0.021D0)
 748       call reada(weightcard,"CTRISS",ctriss,1.001D0)
 749       call reada(weightcard,"DTRISS",dtriss,1.001D0)
 750       write (iout,*) "ATRISS=", atriss
 751       write (iout,*) "BTRISS=", btriss
 752       write (iout,*) "CTRISS=", ctriss
 753       write (iout,*) "DTRISS=", dtriss
 754       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 755       do i=1,maxres
 756         dyn_ss_mask(i)=.false.
 757       enddo
 758       do i=1,maxres-1
 759         do j=i+1,maxres
 760           dyn_ssbond_ij(i,j)=1.0d300
 761         enddo
 762       enddo
 763       call reada(weightcard,"HT",Ht,0.0D0)
 764       if (dyn_ss) then
 765         ss_depth=ebr/wsc-0.25*eps(1,1)
 766         Ht=Ht/wsc-0.25*eps(1,1)
 767         akcm=akcm*wstrain/wsc
 768         akth=akth*wstrain/wsc
 769         akct=akct*wstrain/wsc
 770         v1ss=v1ss*wstrain/wsc
 771         v2ss=v2ss*wstrain/wsc
 772         v3ss=v3ss*wstrain/wsc
 773       else
 774         if (wstrain.ne.0.0) then
 775          ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 776         else
 777           ss_depth=0.0
 778         endif
 779       endif
 780
 781       if(me.eq.king.or..not.out1file) then
 782        write (iout,*) "Parameters of the SS-bond potential:"
 783        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 784      & " AKCT",akct
 785        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 786        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
 787        write (iout,*)" HT",Ht
 788        print *,'indpdb=',indpdb,' pdbref=',pdbref
 789       endif
 790       if (indpdb.gt.0 .or. pdbref) then
 791         read(inp,'(a)') pdbfile
 792         if(me.eq.king.or..not.out1file)
 793      &   write (iout,'(2a)') 'PDB data will be read from file ',
 794      &   pdbfile(:ilen(pdbfile))
 795         open(ipdbin,file=pdbfile,status='old',err=33)
 796         goto 34
 797   33    write (iout,'(a)') 'Error opening PDB file.'
 798         stop
 799   34    continue
 800 c        write (iout,*) 'Begin reading pdb data'
 801 c        call flush(iout)
 802         call readpdb
 803 c        write (iout,*) 'Finished reading pdb data'
 804 c        call flush(iout)
 805         if(me.eq.king.or..not.out1file)
 806      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 807      &   ' nstart_sup=',nstart_sup
 808         do i=1,nres
 809           itype_pdb(i)=itype(i)
 810         enddo
 811         close (ipdbin)
 812         nnt=nstart_sup
 813         nct=nstart_sup+nsup-1
 814         call contact(.false.,ncont_ref,icont_ref,co)
 815
 816         if (sideadd) then
 817          if(me.eq.king.or..not.out1file)
 818      &    write(iout,*)'Adding sidechains'
 819          maxsi=1000
 820          do i=2,nres-1
 821           iti=itype(i)
 822           if (iti.ne.10 .and. itype(i).ne.ntyp1) then
 823             nsi=0
 824             fail=.true.
 825             do while (fail.and.nsi.le.maxsi)
 826               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 827               nsi=nsi+1
 828             enddo
 829             if(fail) write(iout,*)'Adding sidechain failed for res ',
 830      &              i,' after ',nsi,' trials'
 831           endif
 832          enddo
 833         endif
 834       endif
 835       if (indpdb.eq.0) then
 836 C Read sequence if not taken from the pdb file.
 837         read (inp,*) nres
 838 c        print *,'nres=',nres
 839         if (iscode.gt.0) then
 840           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 841         else
 842           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 843         endif
 844 C Convert sequence to numeric code
 845         do i=1,nres
 846           itype(i)=rescode(i,sequence(i),iscode)
 847         enddo
 848 C Assign initial virtual bond lengths
 849         do i=2,nres
 850           vbld(i)=vbl
 851           vbld_inv(i)=vblinv
 852         enddo
 853         do i=2,nres-1
 854           vbld(i+nres)=dsc(iabs(itype(i)))
 855           vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
 856 c          write (iout,*) "i",i," itype",itype(i),
 857 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 858         enddo
 859       endif
 860 c      print *,nres
 861 c      print '(20i4)',(itype(i),i=1,nres)
 862       do i=1,nres
 863 #ifdef PROCOR
 864         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 865 #else
 866         if (itype(i).eq.ntyp1) then
 867 #endif
 868           itel(i)=0
 869 #ifdef PROCOR
 870         else if (iabs(itype(i+1)).ne.20) then
 871 #else
 872         else if (iabs(itype(i)).ne.20) then
 873 #endif
 874           itel(i)=1
 875         else
 876           itel(i)=2
 877         endif
 878       enddo
 879       if(me.eq.king.or..not.out1file)then
 880        write (iout,*) "ITEL"
 881        do i=1,nres-1
 882          write (iout,*) i,itype(i),itel(i)
 883        enddo
 884        print *,'Call Read_Bridge.'
 885       endif
 886       call read_bridge
 887 C 8/13/98 Set limits to generating the dihedral angles
 888       do i=1,nres
 889         phibound(1,i)=-pi
 890         phibound(2,i)=pi
 891       enddo
 892       read (inp,*) ndih_constr
 893       if (ndih_constr.gt.0) then
 894 C        read (inp,*) ftors
 895         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
 896      &  i=1,ndih_constr)
 897         if(me.eq.king.or..not.out1file)then
 898          write (iout,*)
 899      &   'There are',ndih_constr,' constraints on phi angles.'
 900          do i=1,ndih_constr
 901           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
 902      &    ftors(i)
 903          enddo
 904         endif
 905         do i=1,ndih_constr
 906           phi0(i)=deg2rad*phi0(i)
 907           drange(i)=deg2rad*drange(i)
 908         enddo
 909 C        if(me.eq.king.or..not.out1file)
 910 C     &   write (iout,*) 'FTORS',ftors
 911         do i=1,ndih_constr
 912           ii = idih_constr(i)
 913           phibound(1,ii) = phi0(i)-drange(i)
 914           phibound(2,ii) = phi0(i)+drange(i)
 915         enddo
 916       endif
 917 C first setting the theta boundaries to 0 to pi
 918 C this mean that there is no energy penalty for any angle occuring this can be applied
 919 C for generate random conformation but is not implemented in this way
 920 C      do i=1,nres
 921 C        thetabound(1,i)=0
 922 C        thetabound(2,i)=pi
 923 C      enddo
 924 C begin reading theta constrains this is quartic constrains allowing to
 925 C have smooth second derivative
 926       if (with_theta_constr) then
 927 C with_theta_constr is keyword allowing for occurance of theta constrains
 928       read (inp,*) ntheta_constr
 929 C ntheta_constr is the number of theta constrains
 930       if (ntheta_constr.gt.0) then
 931 C        read (inp,*) ftors
 932         read (inp,*) (itheta_constr(i),theta_constr0(i),
 933      &  theta_drange(i),for_thet_constr(i),
 934      &  i=1,ntheta_constr)
 935 C the above code reads from 1 to ntheta_constr
 936 C itheta_constr(i) residue i for which is theta_constr
 937 C theta_constr0 the global minimum value
 938 C theta_drange is range for which there is no energy penalty
 939 C for_thet_constr is the force constant for quartic energy penalty
 940 C E=k*x**4
 941         if(me.eq.king.or..not.out1file)then
 942          write (iout,*)
 943      &   'There are',ntheta_constr,' constraints on phi angles.'
 944          do i=1,ntheta_constr
 945           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
 946      &    theta_drange(i),
 947      &    for_thet_constr(i)
 948          enddo
 949         endif
 950         do i=1,ntheta_constr
 951           theta_constr0(i)=deg2rad*theta_constr0(i)
 952           theta_drange(i)=deg2rad*theta_drange(i)
 953         enddo
 954 C        if(me.eq.king.or..not.out1file)
 955 C     &   write (iout,*) 'FTORS',ftors
 956 C        do i=1,ntheta_constr
 957 C          ii = itheta_constr(i)
 958 C          thetabound(1,ii) = phi0(i)-drange(i)
 959 C          thetabound(2,ii) = phi0(i)+drange(i)
 960 C        enddo
 961       endif ! ntheta_constr.gt.0
 962       endif! with_theta_constr
 963 C
 964 C      with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
 965 C      write (iout,*) "with_dihed_constr ",with_dihed_constr
 966       nnt=1
 967 #ifdef MPI
 968       if (me.eq.king) then
 969 #endif
 970        write (iout,'(a)') 'Boundaries in phi angle sampling:'
 971        do i=1,nres
 972          write (iout,'(a3,i5,2f10.1)')
 973      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 974        enddo
 975 #ifdef MP
 976       endif
 977 #endif
 978       nct=nres
 979 cd      print *,'NNT=',NNT,' NCT=',NCT
 980       if (itype(1).eq.ntyp1) nnt=2
 981       if (itype(nres).eq.ntyp1) nct=nct-1
 982       if (pdbref) then
 983         if(me.eq.king.or..not.out1file)
 984      &   write (iout,'(a,i3)') 'nsup=',nsup
 985         nstart_seq=nnt
 986         if (nsup.le.(nct-nnt+1)) then
 987           do i=0,nct-nnt+1-nsup
 988             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 989               nstart_seq=nnt+i
 990               goto 111
 991             endif
 992           enddo
 993           write (iout,'(a)')
 994      &            'Error - sequences to be superposed do not match.'
 995           stop
 996         else
 997           do i=0,nsup-(nct-nnt+1)
 998             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 999      &      then
1000               nstart_sup=nstart_sup+i
1001               nsup=nct-nnt+1
1002               goto 111
1003             endif
1004           enddo
1005           write (iout,'(a)')
1006      &            'Error - sequences to be superposed do not match.'
1007         endif
1008   111   continue
1009         if (nsup.eq.0) nsup=nct-nnt
1010         if (nstart_sup.eq.0) nstart_sup=nnt
1011         if (nstart_seq.eq.0) nstart_seq=nnt
1012         if(me.eq.king.or..not.out1file)
1013      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1014      &                 ' nstart_seq=',nstart_seq
1015       endif
1016 c--- Zscore rms -------
1017       if (nz_start.eq.0) nz_start=nnt
1018       if (nz_end.eq.0 .and. nsup.gt.0) then
1019         nz_end=nnt+nsup-1
1020       else if (nz_end.eq.0) then
1021         nz_end=nct
1022       endif
1023       if(me.eq.king.or..not.out1file)then
1024        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1025        write (iout,*) 'IZ_SC=',iz_sc
1026       endif
1027 c----------------------
1028       call init_int_table
1029       if (refstr) then
1030         if (.not.pdbref) then
1031           call read_angles(inp,*38)
1032           goto 39
1033    38     write (iout,'(a)') 'Error reading reference structure.'
1034 #ifdef MPI
1035           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1036           stop 'Error reading reference structure'
1037 #endif
1038    39     call chainbuild
1039           call setup_var
1040 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1041           nstart_sup=nnt
1042           nstart_seq=nnt
1043           nsup=nct-nnt+1
1044           kkk=1
1045           do i=1,2*nres
1046             do j=1,3
1047               cref(j,i,kkk)=c(j,i)
1048             enddo
1049           enddo
1050           call contact(.true.,ncont_ref,icont_ref,co)
1051         endif
1052         endif
1053         print *, "A TU"
1054         write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1055         call flush(iout)
1056         if (constr_dist.gt.0) call read_dist_constr
1057         write (iout,*) "After read_dist_constr nhpb",nhpb
1058         if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1059         call hpb_partition
1060         if(me.eq.king.or..not.out1file)
1061      &   write (iout,*) 'Contact order:',co
1062         if (pdbref) then
1063         if(me.eq.king.or..not.out1file)
1064      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1065         do i=1,ncont_ref
1066           do j=1,2
1067             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1068           enddo
1069           if(me.eq.king.or..not.out1file)
1070      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1071      &     icont_ref(1,i),' ',
1072      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1073         enddo
1074         endif
1075 C      endif
1076       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1077      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1078      &    modecalc.ne.10) then
1079 C If input structure hasn't been supplied from the PDB file read or generate
1080 C initial geometry.
1081         if (iranconf.eq.0 .and. .not. extconf) then
1082           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1083      &     write (iout,'(a)') 'Initial geometry will be read in.'
1084           if (read_cart) then
1085             read(inp,'(8f10.5)',end=36,err=36)
1086      &       ((c(l,k),l=1,3),k=1,nres),
1087      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1088             write (iout,*) "Exit READ_CART"
1089             write (iout,'(8f10.5)')
1090      &       ((c(l,k),l=1,3),k=1,nres),
1091      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1092             call int_from_cart1(.true.)
1093             write (iout,*) "Finish INT_TO_CART"
1094             do i=1,nres-1
1095               do j=1,3
1096                 dc(j,i)=c(j,i+1)-c(j,i)
1097                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1098               enddo
1099             enddo
1100             do i=nnt,nct
1101               if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1102                 do j=1,3
1103                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1104                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1105                 enddo
1106               endif
1107             enddo
1108             return
1109           else
1110             call read_angles(inp,*36)
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         else
1136           if(me.eq.king.or..not.out1file)
1137      &     write (iout,'(a)') 'Random-generated initial geometry.'
1138
1139
1140 #ifdef MPI
1141           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1142      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1143 #endif
1144             do itrial=1,100
1145               itmp=1
1146               call gen_rand_conf(itmp,*30)
1147               goto 40
1148    30         write (iout,*) 'Failed to generate random conformation',
1149      &          ', itrial=',itrial
1150               write (*,*) 'Processor:',me,
1151      &          ' Failed to generate random conformation',
1152      &          ' itrial=',itrial
1153               call intout
1154
1155 #ifdef AIX
1156               call flush_(iout)
1157 #else
1158               call flush(iout)
1159 #endif
1160             enddo
1161             write (iout,'(a,i3,a)') 'Processor:',me,
1162      &        ' error in generating random conformation.'
1163             write (*,'(a,i3,a)') 'Processor:',me,
1164      &        ' error in generating random conformation.'
1165             call flush(iout)
1166 #ifdef MPI
1167             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1168    40       continue
1169           endif
1170 #else
1171           write (*,'(a)')
1172      &      ' error in generating random conformation.'
1173           stop
1174    40     continue
1175 #endif
1176         endif
1177       elseif (modecalc.eq.4) then
1178         read (inp,'(a)') intinname
1179         open (intin,file=intinname,status='old',err=333)
1180         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1181      &  write (iout,'(a)') 'intinname',intinname
1182         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1183         goto 334
1184   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1185 #ifdef MPI
1186         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1187 #endif
1188         stop 'Error opening angle file.'
1189   334   continue
1190
1191       endif
1192 C Generate distance constraints, if the PDB structure is to be regularized.
1193       if (nthread.gt.0) then
1194         call read_threadbase
1195       endif
1196       call setup_var
1197       if (me.eq.king .or. .not. out1file)
1198      & call intout
1199       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1200         write (iout,'(/a,i3,a)')
1201      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1202         write (iout,'(20i4)') (iss(i),i=1,ns)
1203        if (dyn_ss) then
1204           write(iout,*)"Running with dynamic disulfide-bond formation"
1205        else
1206         write (iout,'(/a/)') 'Pre-formed links are:'
1207         do i=1,nss
1208           i1=ihpb(i)-nres
1209           i2=jhpb(i)-nres
1210           it1=itype(i1)
1211           it2=itype(i2)
1212           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1213      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1214      &    ebr,forcon(i)
1215         enddo
1216         write (iout,'(a)')
1217        endif
1218       endif
1219       if (ns.gt.0.and.dyn_ss) then
1220           do i=nss+1,nhpb
1221             ihpb(i-nss)=ihpb(i)
1222             jhpb(i-nss)=jhpb(i)
1223             forcon(i-nss)=forcon(i)
1224             dhpb(i-nss)=dhpb(i)
1225           enddo
1226           nhpb=nhpb-nss
1227           nss=0
1228           call hpb_partition
1229           do i=1,ns
1230             dyn_ss_mask(iss(i))=.true.
1231           enddo
1232       endif
1233       if (i2ndstr.gt.0) call secstrp2dihc
1234 c      call geom_to_var(nvar,x)
1235 c      call etotal(energia(0))
1236 c      call enerprint(energia(0))
1237 c      call briefout(0,etot)
1238 c      stop
1239 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1240 cd    write (iout,'(a)') 'Variable list:'
1241 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1242 #ifdef MPI
1243       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1244      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1245      &  'Processor',myrank,': end reading molecular data.'
1246 #endif
1247       print *,"A TU?"
1248       return
1249       end
1250 c--------------------------------------------------------------------------
1251       logical function seq_comp(itypea,itypeb,length)
1252       implicit none
1253       integer length,itypea(length),itypeb(length)
1254       integer i
1255       do i=1,length
1256         if (itypea(i).ne.itypeb(i)) then
1257           seq_comp=.false.
1258           return
1259         endif
1260       enddo
1261       seq_comp=.true.
1262       return
1263       end
1264 c-----------------------------------------------------------------------------
1265       subroutine read_bridge
1266 C Read information about disulfide bridges.
1267       implicit real*8 (a-h,o-z)
1268       include 'DIMENSIONS'
1269 #ifdef MPI
1270       include 'mpif.h'
1271 #endif
1272       include 'COMMON.IOUNITS'
1273       include 'COMMON.GEO'
1274       include 'COMMON.VAR'
1275       include 'COMMON.INTERACT'
1276       include 'COMMON.LOCAL'
1277       include 'COMMON.NAMES'
1278       include 'COMMON.CHAIN'
1279       include 'COMMON.FFIELD'
1280       include 'COMMON.SBRIDGE'
1281       include 'COMMON.HEADER'
1282       include 'COMMON.CONTROL'
1283       include 'COMMON.DBASE'
1284       include 'COMMON.THREAD'
1285       include 'COMMON.TIME1'
1286       include 'COMMON.SETUP'
1287 C Read bridging residues.
1288       read (inp,*) ns,(iss(i),i=1,ns)
1289       print *,'ns=',ns
1290       if(me.eq.king.or..not.out1file)
1291      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1292 C Check whether the specified bridging residues are cystines.
1293       do i=1,ns
1294         if (itype(iss(i)).ne.1) then
1295           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1296      &   'Do you REALLY think that the residue ',
1297      &    restyp(itype(iss(i))),i,
1298      &   ' can form a disulfide bridge?!!!'
1299           write (*,'(2a,i3,a)')
1300      &   'Do you REALLY think that the residue ',
1301      &    restyp(itype(iss(i))),i,
1302      &   ' can form a disulfide bridge?!!!'
1303 #ifdef MPI
1304          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1305          stop
1306 #endif
1307         endif
1308       enddo
1309 C Read preformed bridges.
1310       if (ns.gt.0) then
1311       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1312       if(fg_rank.eq.0)
1313      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1314       if (nss.gt.0) then
1315         nhpb=nss
1316 C Check if the residues involved in bridges are in the specified list of
1317 C bridging residues.
1318         do i=1,nss
1319           do j=1,i-1
1320             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1321      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1322               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1323      &      ' contains residues present in other pairs.'
1324               write (*,'(a,i3,a)') 'Disulfide pair',i,
1325      &      ' contains residues present in other pairs.'
1326 #ifdef MPI
1327               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1328               stop
1329 #endif
1330             endif
1331           enddo
1332           do j=1,ns
1333             if (ihpb(i).eq.iss(j)) goto 10
1334           enddo
1335           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1336    10     continue
1337           do j=1,ns
1338             if (jhpb(i).eq.iss(j)) goto 20
1339           enddo
1340           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1341    20     continue
1342           dhpb(i)=dbr
1343           forcon(i)=fbr
1344         enddo
1345         do i=1,nss
1346           ihpb(i)=ihpb(i)+nres
1347           jhpb(i)=jhpb(i)+nres
1348         enddo
1349       endif
1350       endif
1351       return
1352       end
1353 c----------------------------------------------------------------------------
1354       subroutine read_x(kanal,*)
1355       implicit real*8 (a-h,o-z)
1356       include 'DIMENSIONS'
1357       include 'COMMON.GEO'
1358       include 'COMMON.VAR'
1359       include 'COMMON.CHAIN'
1360       include 'COMMON.IOUNITS'
1361       include 'COMMON.CONTROL'
1362       include 'COMMON.LOCAL'
1363       include 'COMMON.INTERACT'
1364 c Read coordinates from input
1365 c
1366       read(kanal,'(8f10.5)',end=10,err=10)
1367      &  ((c(l,k),l=1,3),k=1,nres),
1368      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1369       do j=1,3
1370         c(j,nres+1)=c(j,1)
1371         c(j,2*nres)=c(j,nres)
1372       enddo
1373       call int_from_cart1(.false.)
1374       do i=1,nres-1
1375         do j=1,3
1376           dc(j,i)=c(j,i+1)-c(j,i)
1377           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1378         enddo
1379       enddo
1380       do i=nnt,nct
1381         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1382           do j=1,3
1383             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1384             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1385           enddo
1386         endif
1387       enddo
1388
1389       return
1390    10 return1
1391       end
1392 c----------------------------------------------------------------------------
1393       subroutine read_threadbase
1394       implicit real*8 (a-h,o-z)
1395       include 'DIMENSIONS'
1396       include 'COMMON.IOUNITS'
1397       include 'COMMON.GEO'
1398       include 'COMMON.VAR'
1399       include 'COMMON.INTERACT'
1400       include 'COMMON.LOCAL'
1401       include 'COMMON.NAMES'
1402       include 'COMMON.CHAIN'
1403       include 'COMMON.FFIELD'
1404       include 'COMMON.SBRIDGE'
1405       include 'COMMON.HEADER'
1406       include 'COMMON.CONTROL'
1407       include 'COMMON.DBASE'
1408       include 'COMMON.THREAD'
1409       include 'COMMON.TIME1'
1410 C Read pattern database for threading.
1411       read (icbase,*) nseq
1412       do i=1,nseq
1413         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1414      &   nres_base(2,i),nres_base(3,i)
1415         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1416      &   nres_base(1,i))
1417 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1418 c    &   nres_base(2,i),nres_base(3,i)
1419 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1420 c    &   nres_base(1,i))
1421       enddo
1422       close (icbase)
1423       if (weidis.eq.0.0D0) weidis=0.1D0
1424       do i=nnt,nct
1425         do j=i+2,nct
1426           nhpb=nhpb+1
1427           ihpb(nhpb)=i
1428           jhpb(nhpb)=j
1429           forcon(nhpb)=weidis
1430         enddo
1431       enddo
1432       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1433       write (iout,'(a,i5)') 'nexcl: ',nexcl
1434       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1435       return
1436       end
1437 c------------------------------------------------------------------------------
1438       subroutine setup_var
1439       implicit real*8 (a-h,o-z)
1440       include 'DIMENSIONS'
1441       include 'COMMON.IOUNITS'
1442       include 'COMMON.GEO'
1443       include 'COMMON.VAR'
1444       include 'COMMON.INTERACT'
1445       include 'COMMON.LOCAL'
1446       include 'COMMON.NAMES'
1447       include 'COMMON.CHAIN'
1448       include 'COMMON.FFIELD'
1449       include 'COMMON.SBRIDGE'
1450       include 'COMMON.HEADER'
1451       include 'COMMON.CONTROL'
1452       include 'COMMON.DBASE'
1453       include 'COMMON.THREAD'
1454       include 'COMMON.TIME1'
1455 C Set up variable list.
1456       ntheta=nres-2
1457       nphi=nres-3
1458       nvar=ntheta+nphi
1459       nside=0
1460       do i=2,nres-1
1461         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1462           nside=nside+1
1463           ialph(i,1)=nvar+nside
1464           ialph(nside,2)=i
1465         endif
1466       enddo
1467       if (indphi.gt.0) then
1468         nvar=nphi
1469       else if (indback.gt.0) then
1470         nvar=nphi+ntheta
1471       else
1472         nvar=nvar+2*nside
1473       endif
1474 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1475       return
1476       end
1477 c----------------------------------------------------------------------------
1478       subroutine gen_dist_constr
1479 C Generate CA distance constraints.
1480       implicit real*8 (a-h,o-z)
1481       include 'DIMENSIONS'
1482       include 'COMMON.IOUNITS'
1483       include 'COMMON.GEO'
1484       include 'COMMON.VAR'
1485       include 'COMMON.INTERACT'
1486       include 'COMMON.LOCAL'
1487       include 'COMMON.NAMES'
1488       include 'COMMON.CHAIN'
1489       include 'COMMON.FFIELD'
1490       include 'COMMON.SBRIDGE'
1491       include 'COMMON.HEADER'
1492       include 'COMMON.CONTROL'
1493       include 'COMMON.DBASE'
1494       include 'COMMON.THREAD'
1495       include 'COMMON.TIME1'
1496       dimension itype_pdb(maxres)
1497       common /pizda/ itype_pdb
1498       character*2 iden
1499 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1500 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1501 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1502 cd     & ' nsup',nsup
1503       do i=nstart_sup,nstart_sup+nsup-1
1504 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1505 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1506         do j=i+2,nstart_sup+nsup-1
1507           nhpb=nhpb+1
1508           ihpb(nhpb)=i+nstart_seq-nstart_sup
1509           jhpb(nhpb)=j+nstart_seq-nstart_sup
1510           forcon(nhpb)=weidis
1511           dhpb(nhpb)=dist(i,j)
1512         enddo
1513       enddo
1514 cd      write (iout,'(a)') 'Distance constraints:'
1515 cd      do i=nss+1,nhpb
1516 cd        ii=ihpb(i)
1517 cd        jj=jhpb(i)
1518 cd        iden='CA'
1519 cd        if (ii.gt.nres) then
1520 cd          iden='SC'
1521 cd          ii=ii-nres
1522 cd          jj=jj-nres
1523 cd        endif
1524 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1525 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1526 cd     &  dhpb(i),forcon(i)
1527 cd      enddo
1528       return
1529       end
1530 c----------------------------------------------------------------------------
1531       subroutine map_read
1532       implicit real*8 (a-h,o-z)
1533       include 'DIMENSIONS'
1534       include 'COMMON.MAP'
1535       include 'COMMON.IOUNITS'
1536       character*3 angid(4) /'THE','PHI','ALP','OME'/
1537       character*80 mapcard,ucase
1538       do imap=1,nmap
1539         read (inp,'(a)') mapcard
1540         mapcard=ucase(mapcard)
1541         if (index(mapcard,'PHI').gt.0) then
1542           kang(imap)=1
1543         else if (index(mapcard,'THE').gt.0) then
1544           kang(imap)=2
1545         else if (index(mapcard,'ALP').gt.0) then
1546           kang(imap)=3
1547         else if (index(mapcard,'OME').gt.0) then
1548           kang(imap)=4
1549         else
1550           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1551           stop 'Error - illegal variable spec in MAP card.'
1552         endif
1553         call readi (mapcard,'RES1',res1(imap),0)
1554         call readi (mapcard,'RES2',res2(imap),0)
1555         if (res1(imap).eq.0) then
1556           res1(imap)=res2(imap)
1557         else if (res2(imap).eq.0) then
1558           res2(imap)=res1(imap)
1559         endif
1560         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1561           write (iout,'(a)')
1562      &    'Error - illegal definition of variable group in MAP.'
1563           stop 'Error - illegal definition of variable group in MAP.'
1564         endif
1565         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1566         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1567         call readi(mapcard,'NSTEP',nstep(imap),0)
1568         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1569           write (iout,'(a)')
1570      &     'Illegal boundary and/or step size specification in MAP.'
1571           stop 'Illegal boundary and/or step size specification in MAP.'
1572         endif
1573       enddo ! imap
1574       return
1575       end
1576 c----------------------------------------------------------------------------
1577       subroutine csaread
1578       implicit real*8 (a-h,o-z)
1579       include 'DIMENSIONS'
1580       include 'COMMON.IOUNITS'
1581       include 'COMMON.GEO'
1582       include 'COMMON.CSA'
1583       include 'COMMON.BANK'
1584       include 'COMMON.CONTROL'
1585       character*80 ucase
1586       character*620 mcmcard
1587       call card_concat(mcmcard)
1588
1589       call readi(mcmcard,'NCONF',nconf,50)
1590       call readi(mcmcard,'NADD',nadd,0)
1591       call readi(mcmcard,'JSTART',jstart,1)
1592       call readi(mcmcard,'JEND',jend,1)
1593       call readi(mcmcard,'NSTMAX',nstmax,500000)
1594       call readi(mcmcard,'N0',n0,1)
1595       call readi(mcmcard,'N1',n1,6)
1596       call readi(mcmcard,'N2',n2,4)
1597       call readi(mcmcard,'N3',n3,0)
1598       call readi(mcmcard,'N4',n4,0)
1599       call readi(mcmcard,'N5',n5,0)
1600       call readi(mcmcard,'N6',n6,10)
1601       call readi(mcmcard,'N7',n7,0)
1602       call readi(mcmcard,'N8',n8,0)
1603       call readi(mcmcard,'N9',n9,0)
1604       call readi(mcmcard,'N14',n14,0)
1605       call readi(mcmcard,'N15',n15,0)
1606       call readi(mcmcard,'N16',n16,0)
1607       call readi(mcmcard,'N17',n17,0)
1608       call readi(mcmcard,'N18',n18,0)
1609
1610       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1611
1612       call readi(mcmcard,'NDIFF',ndiff,2)
1613       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1614       call readi(mcmcard,'IS1',is1,1)
1615       call readi(mcmcard,'IS2',is2,8)
1616       call readi(mcmcard,'NRAN0',nran0,4)
1617       call readi(mcmcard,'NRAN1',nran1,2)
1618       call readi(mcmcard,'IRR',irr,1)
1619       call readi(mcmcard,'NSEED',nseed,20)
1620       call readi(mcmcard,'NTOTAL',ntotal,10000)
1621       call reada(mcmcard,'CUT1',cut1,2.0d0)
1622       call reada(mcmcard,'CUT2',cut2,5.0d0)
1623       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1624       call readi(mcmcard,'ICMAX',icmax,3)
1625       call readi(mcmcard,'IRESTART',irestart,0)
1626 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1627       ntbankm=0
1628 c!bankt
1629       call reada(mcmcard,'DELE',dele,20.0d0)
1630       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1631       call readi(mcmcard,'IREF',iref,0)
1632       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1633       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1634       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1635       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1636       write (iout,*) "NCONF_IN",nconf_in
1637       return
1638       end
1639 c----------------------------------------------------------------------------
1640 cfmc      subroutine mcmfread
1641 cfmc      implicit real*8 (a-h,o-z)
1642 cfmc      include 'DIMENSIONS'
1643 cfmc      include 'COMMON.MCMF'
1644 cfmc      include 'COMMON.IOUNITS'
1645 cfmc      include 'COMMON.GEO'
1646 cfmc      character*80 ucase
1647 cfmc      character*620 mcmcard
1648 cfmc      call card_concat(mcmcard)
1649 cfmc
1650 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1651 cfmc      write(iout,*)'MAXRANT=',maxrant
1652 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1653 cfmc      write(iout,*)'MAXFAM=',maxfam
1654 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1655 cfmc      write(iout,*)'NNET1=',nnet1
1656 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1657 cfmc      write(iout,*)'NNET2=',nnet2
1658 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1659 cfmc      write(iout,*)'NNET3=',nnet3
1660 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1661 cfmc      write(iout,*)'ILASTT=',ilastt
1662 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1663 cfmc      write(iout,*)'MAXSTR=',maxstr
1664 cfmc      maxstr_f=maxstr/maxfam
1665 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1666 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1667 cfmc      write(iout,*)'NMCMF=',nmcmf
1668 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1669 cfmc      write(iout,*)'IFOCUS=',ifocus
1670 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1671 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1672 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1673 cfmc      write(iout,*)'INTPRT=',intprt
1674 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1675 cfmc      write(iout,*)'IPRT=',iprt
1676 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1677 cfmc      write(iout,*)'IMAXTR=',imaxtr
1678 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1679 cfmc      write(iout,*)'MAXEVEN=',maxeven
1680 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1681 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1682 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1683 cfmc      write(iout,*)'INIMIN=',inimin
1684 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1685 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1686 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1687 cfmc      write(iout,*)'NTHREAD=',nthread
1688 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1689 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1690 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1691 cfmc      write(iout,*)'MAXPERT=',maxpert
1692 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1693 cfmc      write(iout,*)'IRMSD=',irmsd
1694 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1695 cfmc      write(iout,*)'DENEMIN=',denemin
1696 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1697 cfmc      write(iout,*)'RCUT1S=',rcut1s
1698 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1699 cfmc      write(iout,*)'RCUT1E=',rcut1e
1700 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1701 cfmc      write(iout,*)'RCUT2S=',rcut2s
1702 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1703 cfmc      write(iout,*)'RCUT2E=',rcut2e
1704 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1705 cfmc      write(iout,*)'DPERT1=',d_pert1
1706 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1707 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1708 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1709 cfmc      write(iout,*)'DPERT2=',d_pert2
1710 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1711 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1712 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1713 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1714 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1715 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1716 cfmc      d_pert1=deg2rad*d_pert1
1717 cfmc      d_pert1a=deg2rad*d_pert1a
1718 cfmc      d_pert2=deg2rad*d_pert2
1719 cfmc      d_pert2a=deg2rad*d_pert2a
1720 cfmc      d_pert2b=deg2rad*d_pert2b
1721 cfmc      d_pert2c=deg2rad*d_pert2c
1722 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1723 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1724 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1725 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1726 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1727 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1728 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1729 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1730 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1731 cfmc      write(iout,*)'RCUTINI=',rcutini
1732 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1733 cfmc      write(iout,*)'GRAT=',grat
1734 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1735 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1736 cfmc
1737 cfmc      return
1738 cfmc      end
1739 c----------------------------------------------------------------------------
1740       subroutine mcmread
1741       implicit real*8 (a-h,o-z)
1742       include 'DIMENSIONS'
1743       include 'COMMON.MCM'
1744       include 'COMMON.MCE'
1745       include 'COMMON.IOUNITS'
1746       character*80 ucase
1747       character*320 mcmcard
1748       call card_concat(mcmcard)
1749       call readi(mcmcard,'MAXACC',maxacc,100)
1750       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1751       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1752       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1753       call readi(mcmcard,'MAXREPM',maxrepm,200)
1754       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1755       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1756       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1757       call reada(mcmcard,'E_UP',e_up,5.0D0)
1758       call reada(mcmcard,'DELTE',delte,0.1D0)
1759       call readi(mcmcard,'NSWEEP',nsweep,5)
1760       call readi(mcmcard,'NSTEPH',nsteph,0)
1761       call readi(mcmcard,'NSTEPC',nstepc,0)
1762       call reada(mcmcard,'TMIN',tmin,298.0D0)
1763       call reada(mcmcard,'TMAX',tmax,298.0D0)
1764       call readi(mcmcard,'NWINDOW',nwindow,0)
1765       call readi(mcmcard,'PRINT_MC',print_mc,0)
1766       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1767       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1768       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1769       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1770       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1771       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1772       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1773       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1774       if (nwindow.gt.0) then
1775         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1776         do i=1,nwindow
1777           winlen(i)=winend(i)-winstart(i)+1
1778         enddo
1779       endif
1780       if (tmax.lt.tmin) tmax=tmin
1781       if (tmax.eq.tmin) then
1782         nstepc=0
1783         nsteph=0
1784       endif
1785       if (nstepc.gt.0 .and. nsteph.gt.0) then
1786         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1787         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1788       endif
1789 C Probabilities of different move types
1790       sumpro_type(0)=0.0D0
1791       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1792       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1793       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1794       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1795       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1796       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1797       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1798       do i=1,MaxMoveType
1799         print *,'i',i,' sumprotype',sumpro_type(i)
1800         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1801         print *,'i',i,' sumprotype',sumpro_type(i)
1802       enddo
1803       return
1804       end
1805 c----------------------------------------------------------------------------
1806       subroutine read_minim
1807       implicit real*8 (a-h,o-z)
1808       include 'DIMENSIONS'
1809       include 'COMMON.MINIM'
1810       include 'COMMON.IOUNITS'
1811       character*80 ucase
1812       character*320 minimcard
1813       call card_concat(minimcard)
1814       call readi(minimcard,'MAXMIN',maxmin,2000)
1815       call readi(minimcard,'MAXFUN',maxfun,5000)
1816       call readi(minimcard,'MINMIN',minmin,maxmin)
1817       call readi(minimcard,'MINFUN',minfun,maxmin)
1818       call reada(minimcard,'TOLF',tolf,1.0D-2)
1819       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1820       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1821       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1822       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1823       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1824      &         'Options in energy minimization:'
1825       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1826      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1827      & 'MinMin:',MinMin,' MinFun:',MinFun,
1828      & ' TolF:',TolF,' RTolF:',RTolF
1829       return
1830       end
1831 c----------------------------------------------------------------------------
1832       subroutine read_angles(kanal,*)
1833       implicit real*8 (a-h,o-z)
1834       include 'DIMENSIONS'
1835       include 'COMMON.GEO'
1836       include 'COMMON.VAR'
1837       include 'COMMON.CHAIN'
1838       include 'COMMON.IOUNITS'
1839       include 'COMMON.CONTROL'
1840 c Read angles from input
1841 c
1842        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1843        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1844        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1845        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1846
1847        do i=1,nres
1848 c 9/7/01 avoid 180 deg valence angle
1849         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1850 c
1851         theta(i)=deg2rad*theta(i)
1852         phi(i)=deg2rad*phi(i)
1853         alph(i)=deg2rad*alph(i)
1854         omeg(i)=deg2rad*omeg(i)
1855        enddo
1856       return
1857    10 return1
1858       end
1859 c----------------------------------------------------------------------------
1860       subroutine reada(rekord,lancuch,wartosc,default)
1861       implicit none
1862       character*(*) rekord,lancuch
1863       double precision wartosc,default
1864       integer ilen,iread
1865       external ilen
1866       iread=index(rekord,lancuch)
1867       if (iread.eq.0) then
1868         wartosc=default
1869         return
1870       endif
1871       iread=iread+ilen(lancuch)+1
1872       read (rekord(iread:),*,err=10,end=10) wartosc
1873       return
1874   10  wartosc=default
1875       return
1876       end
1877 c----------------------------------------------------------------------------
1878       subroutine readi(rekord,lancuch,wartosc,default)
1879       implicit none
1880       character*(*) rekord,lancuch
1881       integer wartosc,default
1882       integer ilen,iread
1883       external ilen
1884       iread=index(rekord,lancuch)
1885       if (iread.eq.0) then
1886         wartosc=default
1887         return
1888       endif
1889       iread=iread+ilen(lancuch)+1
1890       read (rekord(iread:),*,err=10,end=10) wartosc
1891       return
1892   10  wartosc=default
1893       return
1894       end
1895 c----------------------------------------------------------------------------
1896       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1897       implicit none
1898       integer dim,i
1899       integer tablica(dim),default
1900       character*(*) rekord,lancuch
1901       character*80 aux
1902       integer ilen,iread
1903       external ilen
1904       do i=1,dim
1905         tablica(i)=default
1906       enddo
1907       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1908       if (iread.eq.0) return
1909       iread=iread+ilen(lancuch)+1
1910       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1911    10 return
1912       end
1913 c----------------------------------------------------------------------------
1914       subroutine multreada(rekord,lancuch,tablica,dim,default)
1915       implicit none
1916       integer dim,i
1917       double precision tablica(dim),default
1918       character*(*) rekord,lancuch
1919       character*80 aux
1920       integer ilen,iread
1921       external ilen
1922       do i=1,dim
1923         tablica(i)=default
1924       enddo
1925       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1926       if (iread.eq.0) return
1927       iread=iread+ilen(lancuch)+1
1928       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1929    10 return
1930       end
1931 c----------------------------------------------------------------------------
1932       subroutine openunits
1933       implicit real*8 (a-h,o-z)
1934       include 'DIMENSIONS'
1935 #ifdef MPI
1936       include 'mpif.h'
1937       character*16 form,nodename
1938       integer nodelen
1939 #endif
1940       include 'COMMON.SETUP'
1941       include 'COMMON.IOUNITS'
1942       include 'COMMON.MD'
1943       include 'COMMON.CONTROL'
1944       integer lenpre,lenpot,ilen,lentmp
1945       external ilen
1946       character*3 out1file_text,ucase
1947       character*3 ll
1948       external ucase
1949 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1950       call getenv_loc("PREFIX",prefix)
1951       pref_orig = prefix
1952       call getenv_loc("POT",pot)
1953       call getenv_loc("DIRTMP",tmpdir)
1954       call getenv_loc("CURDIR",curdir)
1955       call getenv_loc("OUT1FILE",out1file_text)
1956 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1957       out1file_text=ucase(out1file_text)
1958       if (out1file_text(1:1).eq."Y") then
1959         out1file=.true.
1960       else
1961         out1file=fg_rank.gt.0
1962       endif
1963       lenpre=ilen(prefix)
1964       lenpot=ilen(pot)
1965       lentmp=ilen(tmpdir)
1966       if (lentmp.gt.0) then
1967           write (*,'(80(1h!))')
1968           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1969           write (*,'(80(1h!))')
1970           write (*,*)"All output files will be on node /tmp directory."
1971 #ifdef MPI
1972         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1973         if (me.eq.king) then
1974           write (*,*) "The master node is ",nodename
1975         else if (fg_rank.eq.0) then
1976           write (*,*) "I am the CG slave node ",nodename
1977         else
1978           write (*,*) "I am the FG slave node ",nodename
1979         endif
1980 #endif
1981         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1982         lenpre = lentmp+lenpre+1
1983       endif
1984       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1985 C Get the names and open the input files
1986 #if defined(WINIFL) || defined(WINPGI)
1987       open(1,file=pref_orig(:ilen(pref_orig))//
1988      &  '.inp',status='old',readonly,shared)
1989        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1990 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1991 C Get parameter filenames and open the parameter files.
1992       call getenv_loc('BONDPAR',bondname)
1993       open (ibond,file=bondname,status='old',readonly,shared)
1994       call getenv_loc('THETPAR',thetname)
1995       open (ithep,file=thetname,status='old',readonly,shared)
1996       call getenv_loc('ROTPAR',rotname)
1997       open (irotam,file=rotname,status='old',readonly,shared)
1998       call getenv_loc('TORPAR',torname)
1999       open (itorp,file=torname,status='old',readonly,shared)
2000       call getenv_loc('TORDPAR',tordname)
2001       open (itordp,file=tordname,status='old',readonly,shared)
2002       call getenv_loc('FOURIER',fouriername)
2003       open (ifourier,file=fouriername,status='old',readonly,shared)
2004       call getenv_loc('ELEPAR',elename)
2005       open (ielep,file=elename,status='old',readonly,shared)
2006       call getenv_loc('SIDEPAR',sidename)
2007       open (isidep,file=sidename,status='old',readonly,shared)
2008 #elif (defined CRAY) || (defined AIX)
2009       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2010      &  action='read')
2011 c      print *,"Processor",myrank," opened file 1"
2012       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2013 c      print *,"Processor",myrank," opened file 9"
2014 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2015 C Get parameter filenames and open the parameter files.
2016       call getenv_loc('BONDPAR',bondname)
2017       open (ibond,file=bondname,status='old',action='read')
2018 c      print *,"Processor",myrank," opened file IBOND"
2019       call getenv_loc('THETPAR',thetname)
2020       open (ithep,file=thetname,status='old',action='read')
2021 c      print *,"Processor",myrank," opened file ITHEP"
2022       call getenv_loc('ROTPAR',rotname)
2023       open (irotam,file=rotname,status='old',action='read')
2024 c      print *,"Processor",myrank," opened file IROTAM"
2025       call getenv_loc('TORPAR',torname)
2026       open (itorp,file=torname,status='old',action='read')
2027 c      print *,"Processor",myrank," opened file ITORP"
2028       call getenv_loc('TORDPAR',tordname)
2029       open (itordp,file=tordname,status='old',action='read')
2030 c      print *,"Processor",myrank," opened file ITORDP"
2031       call getenv_loc('SCCORPAR',sccorname)
2032       open (isccor,file=sccorname,status='old',action='read')
2033 c      print *,"Processor",myrank," opened file ISCCOR"
2034       call getenv_loc('FOURIER',fouriername)
2035       open (ifourier,file=fouriername,status='old',action='read')
2036 c      print *,"Processor",myrank," opened file IFOURIER"
2037       call getenv_loc('ELEPAR',elename)
2038       open (ielep,file=elename,status='old',action='read')
2039 c      print *,"Processor",myrank," opened file IELEP"
2040       call getenv_loc('SIDEPAR',sidename)
2041       open (isidep,file=sidename,status='old',action='read')
2042 c      print *,"Processor",myrank," opened file ISIDEP"
2043 c      print *,"Processor",myrank," opened parameter files"
2044 #elif (defined G77)
2045       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2046       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2047 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2048 C Get parameter filenames and open the parameter files.
2049       call getenv_loc('BONDPAR',bondname)
2050       open (ibond,file=bondname,status='old')
2051       call getenv_loc('THETPAR',thetname)
2052       open (ithep,file=thetname,status='old')
2053       call getenv_loc('ROTPAR',rotname)
2054       open (irotam,file=rotname,status='old')
2055       call getenv_loc('TORPAR',torname)
2056       open (itorp,file=torname,status='old')
2057       call getenv_loc('TORDPAR',tordname)
2058       open (itordp,file=tordname,status='old')
2059       call getenv_loc('SCCORPAR',sccorname)
2060       open (isccor,file=sccorname,status='old')
2061       call getenv_loc('FOURIER',fouriername)
2062       open (ifourier,file=fouriername,status='old')
2063       call getenv_loc('ELEPAR',elename)
2064       open (ielep,file=elename,status='old')
2065       call getenv_loc('SIDEPAR',sidename)
2066       open (isidep,file=sidename,status='old')
2067 #else
2068       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2069      &  readonly)
2070        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2071 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2072 C Get parameter filenames and open the parameter files.
2073       call getenv_loc('BONDPAR',bondname)
2074       open (ibond,file=bondname,status='old',readonly)
2075       call getenv_loc('THETPAR',thetname)
2076       open (ithep,file=thetname,status='old',readonly)
2077       call getenv_loc('ROTPAR',rotname)
2078       open (irotam,file=rotname,status='old',readonly)
2079       call getenv_loc('TORPAR',torname)
2080       open (itorp,file=torname,status='old',readonly)
2081       call getenv_loc('TORDPAR',tordname)
2082       open (itordp,file=tordname,status='old',readonly)
2083       call getenv_loc('SCCORPAR',sccorname)
2084       open (isccor,file=sccorname,status='old',readonly)
2085 #ifndef CRYST_THETA
2086       call getenv_loc('THETPARPDB',thetname_pdb)
2087       print *,"thetname_pdb ",thetname_pdb
2088       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2089       print *,ithep_pdb," opened"
2090 #endif
2091       call getenv_loc('FOURIER',fouriername)
2092       open (ifourier,file=fouriername,status='old',readonly)
2093       call getenv_loc('ELEPAR',elename)
2094       open (ielep,file=elename,status='old',readonly)
2095       call getenv_loc('SIDEPAR',sidename)
2096       open (isidep,file=sidename,status='old',readonly)
2097       call getenv_loc('LIPTRANPAR',liptranname)
2098       open (iliptranpar,file=liptranname,status='old',action='read')
2099 #ifndef CRYST_SC
2100       call getenv_loc('ROTPARPDB',rotname_pdb)
2101       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2102 #endif
2103 #endif
2104 #ifndef OLDSCP
2105 C
2106 C 8/9/01 In the newest version SCp interaction constants are read from a file
2107 C Use -DOLDSCP to use hard-coded constants instead.
2108 C
2109       call getenv_loc('SCPPAR',scpname)
2110 #if defined(WINIFL) || defined(WINPGI)
2111       open (iscpp,file=scpname,status='old',readonly,shared)
2112 #elif (defined CRAY)  || (defined AIX)
2113       open (iscpp,file=scpname,status='old',action='read')
2114 #elif (defined G77)
2115       open (iscpp,file=scpname,status='old')
2116 #else
2117       open (iscpp,file=scpname,status='old',readonly)
2118 #endif
2119 #endif
2120       call getenv_loc('PATTERN',patname)
2121 #if defined(WINIFL) || defined(WINPGI)
2122       open (icbase,file=patname,status='old',readonly,shared)
2123 #elif (defined CRAY)  || (defined AIX)
2124       open (icbase,file=patname,status='old',action='read')
2125 #elif (defined G77)
2126       open (icbase,file=patname,status='old')
2127 #else
2128       open (icbase,file=patname,status='old',readonly)
2129 #endif
2130 #ifdef MPI
2131 C Open output file only for CG processes
2132 c      print *,"Processor",myrank," fg_rank",fg_rank
2133       if (fg_rank.eq.0) then
2134
2135       if (nodes.eq.1) then
2136         npos=3
2137       else
2138         npos = dlog10(dfloat(nodes-1))+1
2139       endif
2140       if (npos.lt.3) npos=3
2141       write (liczba,'(i1)') npos
2142       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2143      &  //')'
2144       write (liczba,form) me
2145       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2146      &  liczba(:ilen(liczba))
2147       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2148      &  //'.int'
2149       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2150      &  //'.pdb'
2151       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2152      &  liczba(:ilen(liczba))//'.mol2'
2153       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2154      &  liczba(:ilen(liczba))//'.stat'
2155       if (lentmp.gt.0)
2156      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2157      &      //liczba(:ilen(liczba))//'.stat')
2158       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2159      &  //'.rst'
2160       if(usampl) then
2161           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2162      & liczba(:ilen(liczba))//'.const'
2163       endif
2164
2165       endif
2166 #else
2167       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2168       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2169       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2170       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2171       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2172       if (lentmp.gt.0)
2173      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2174      &    //'.stat')
2175       rest2name=prefix(:ilen(prefix))//'.rst'
2176       if(usampl) then
2177          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2178       endif
2179 #endif
2180 #if defined(AIX) || defined(PGI)
2181       if (me.eq.king .or. .not. out1file)
2182      &   open(iout,file=outname,status='unknown')
2183 #ifdef DEBUG
2184       if (fg_rank.gt.0) then
2185         write (liczba,'(i3.3)') myrank/nfgtasks
2186         write (ll,'(bz,i3.3)') fg_rank
2187         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2188      &   status='unknown')
2189       endif
2190 #endif
2191       if(me.eq.king) then
2192        open(igeom,file=intname,status='unknown',position='append')
2193        open(ipdb,file=pdbname,status='unknown')
2194        open(imol2,file=mol2name,status='unknown')
2195        open(istat,file=statname,status='unknown',position='append')
2196       else
2197 c1out       open(iout,file=outname,status='unknown')
2198       endif
2199 #else
2200       if (me.eq.king .or. .not.out1file)
2201      &    open(iout,file=outname,status='unknown')
2202 #ifdef DEBUG
2203       if (fg_rank.gt.0) then
2204         write (liczba,'(i3.3)') myrank/nfgtasks
2205         write (ll,'(bz,i3.3)') fg_rank
2206         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2207      &   status='unknown')
2208       endif
2209 #endif
2210       if(me.eq.king) then
2211        open(igeom,file=intname,status='unknown',access='append')
2212        open(ipdb,file=pdbname,status='unknown')
2213        open(imol2,file=mol2name,status='unknown')
2214        open(istat,file=statname,status='unknown',access='append')
2215       else
2216 c1out       open(iout,file=outname,status='unknown')
2217       endif
2218 #endif
2219       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2220       csa_seed=prefix(:lenpre)//'.CSA.seed'
2221       csa_history=prefix(:lenpre)//'.CSA.history'
2222       csa_bank=prefix(:lenpre)//'.CSA.bank'
2223       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2224       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2225       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2226 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2227       csa_int=prefix(:lenpre)//'.int'
2228       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2229       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2230       csa_in=prefix(:lenpre)//'.CSA.in'
2231 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2232 C Write file names
2233       if (me.eq.king)then
2234       write (iout,'(80(1h-))')
2235       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2236       write (iout,'(80(1h-))')
2237       write (iout,*) "Input file                      : ",
2238      &  pref_orig(:ilen(pref_orig))//'.inp'
2239       write (iout,*) "Output file                     : ",
2240      &  outname(:ilen(outname))
2241       write (iout,*)
2242       write (iout,*) "Sidechain potential file        : ",
2243      &  sidename(:ilen(sidename))
2244 #ifndef OLDSCP
2245       write (iout,*) "SCp potential file              : ",
2246      &  scpname(:ilen(scpname))
2247 #endif
2248       write (iout,*) "Electrostatic potential file    : ",
2249      &  elename(:ilen(elename))
2250       write (iout,*) "Cumulant coefficient file       : ",
2251      &  fouriername(:ilen(fouriername))
2252       write (iout,*) "Torsional parameter file        : ",
2253      &  torname(:ilen(torname))
2254       write (iout,*) "Double torsional parameter file : ",
2255      &  tordname(:ilen(tordname))
2256       write (iout,*) "SCCOR parameter file : ",
2257      &  sccorname(:ilen(sccorname))
2258       write (iout,*) "Bond & inertia constant file    : ",
2259      &  bondname(:ilen(bondname))
2260       write (iout,*) "Bending parameter file          : ",
2261      &  thetname(:ilen(thetname))
2262       write (iout,*) "Rotamer parameter file          : ",
2263      &  rotname(:ilen(rotname))
2264       write (iout,*) "Thetpdb parameter file          : ",
2265      &  thetname_pdb(:ilen(thetname_pdb))
2266       write (iout,*) "Threading database              : ",
2267      &  patname(:ilen(patname))
2268       if (lentmp.ne.0)
2269      &write (iout,*)" DIRTMP                          : ",
2270      &  tmpdir(:lentmp)
2271       write (iout,'(80(1h-))')
2272       endif
2273       return
2274       end
2275 c----------------------------------------------------------------------------
2276       subroutine card_concat(card)
2277       implicit real*8 (a-h,o-z)
2278       include 'DIMENSIONS'
2279       include 'COMMON.IOUNITS'
2280       character*(*) card
2281       character*80 karta,ucase
2282       external ilen
2283       read (inp,'(a)') karta
2284       karta=ucase(karta)
2285       card=' '
2286       do while (karta(80:80).eq.'&')
2287         card=card(:ilen(card)+1)//karta(:79)
2288         read (inp,'(a)') karta
2289         karta=ucase(karta)
2290       enddo
2291       card=card(:ilen(card)+1)//karta
2292       return
2293       end
2294 c----------------------------------------------------------------------------------
2295       subroutine readrst
2296       implicit real*8 (a-h,o-z)
2297       include 'DIMENSIONS'
2298       include 'COMMON.CHAIN'
2299       include 'COMMON.IOUNITS'
2300       include 'COMMON.MD'
2301       open(irest2,file=rest2name,status='unknown')
2302       read(irest2,*) totT,EK,potE,totE,t_bath
2303       totTafm=totT
2304       do i=1,2*nres
2305          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2306       enddo
2307       do i=1,2*nres
2308          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2309       enddo
2310       if(usampl) then
2311              read (irest2,*) iset
2312       endif
2313       close(irest2)
2314       return
2315       end
2316 c---------------------------------------------------------------------------------
2317       subroutine read_fragments
2318       implicit real*8 (a-h,o-z)
2319       include 'DIMENSIONS'
2320 #ifdef MPI
2321       include 'mpif.h'
2322 #endif
2323       include 'COMMON.SETUP'
2324       include 'COMMON.CHAIN'
2325       include 'COMMON.IOUNITS'
2326       include 'COMMON.MD'
2327       include 'COMMON.CONTROL'
2328       read(inp,*) nset,nfrag,npair,nfrag_back
2329       if(me.eq.king.or..not.out1file)
2330      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2331      &  " nfrag_back",nfrag_back
2332       do iset=1,nset
2333          read(inp,*) mset(iset)
2334        do i=1,nfrag
2335          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2336      &     qinfrag(i,iset)
2337          if(me.eq.king.or..not.out1file)
2338      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2339      &     ifrag(2,i,iset), qinfrag(i,iset)
2340        enddo
2341        do i=1,npair
2342         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2343      &    qinpair(i,iset)
2344         if(me.eq.king.or..not.out1file)
2345      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2346      &    ipair(2,i,iset), qinpair(i,iset)
2347        enddo
2348        do i=1,nfrag_back
2349         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2350      &     wfrag_back(3,i,iset),
2351      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2352         if(me.eq.king.or..not.out1file)
2353      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2354      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2355        enddo
2356       enddo
2357       return
2358       end
2359 C---------------------------------------------------------------------------
2360       subroutine read_afminp
2361             implicit real*8 (a-h,o-z)
2362       include 'DIMENSIONS'
2363 #ifdef MPI
2364       include 'mpif.h'
2365 #endif
2366       include 'COMMON.SETUP'
2367       include 'COMMON.CONTROL'
2368       include 'COMMON.CHAIN'
2369       include 'COMMON.IOUNITS'
2370       include 'COMMON.SBRIDGE'
2371       character*320 afmcard
2372       print *, "wchodze"
2373       call card_concat(afmcard)
2374       call readi(afmcard,"BEG",afmbeg,0)
2375       call readi(afmcard,"END",afmend,0)
2376       call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2377       call reada(afmcard,"VEL",velAFMconst,0.0d0)
2378       print *,'FORCE=' ,forceAFMconst
2379 CCCC NOW PROPERTIES FOR AFM
2380        distafminit=0.0d0
2381        do i=1,3
2382         distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2383        enddo
2384         distafminit=dsqrt(distafminit)
2385         print *,'initdist',distafminit
2386       return
2387       end
2388
2389 c-------------------------------------------------------------------------------
2390       subroutine read_dist_constr
2391       implicit real*8 (a-h,o-z)
2392       include 'DIMENSIONS'
2393 #ifdef MPI
2394       include 'mpif.h'
2395 #endif
2396       include 'COMMON.SETUP'
2397       include 'COMMON.CONTROL'
2398       include 'COMMON.CHAIN'
2399       include 'COMMON.IOUNITS'
2400       include 'COMMON.SBRIDGE'
2401       integer ifrag_(2,100),ipair_(2,100)
2402       double precision wfrag_(100),wpair_(100)
2403       character*500 controlcard
2404       print *, "WCHODZE"
2405       write (iout,*) "Calling read_dist_constr"
2406 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2407 c      call flush(iout)
2408       call card_concat(controlcard)
2409       call readi(controlcard,"NFRAG",nfrag_,0)
2410       call readi(controlcard,"NPAIR",npair_,0)
2411       call readi(controlcard,"NDIST",ndist_,0)
2412       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2413       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2414       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2415       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2416       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2417 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2418 c      write (iout,*) "IFRAG"
2419 c      do i=1,nfrag_
2420 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2421 c      enddo
2422 c      write (iout,*) "IPAIR"
2423 c      do i=1,npair_
2424 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2425 c      enddo
2426       call flush(iout)
2427       do i=1,nfrag_
2428         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2429         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2430      &    ifrag_(2,i)=nstart_sup+nsup-1
2431 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2432         call flush(iout)
2433         if (wfrag_(i).gt.0.0d0) then
2434         do j=ifrag_(1,i),ifrag_(2,i)-1
2435           do k=j+1,ifrag_(2,i)
2436 c            write (iout,*) "j",j," k",k
2437             ddjk=dist(j,k)
2438             if (constr_dist.eq.1) then
2439             nhpb=nhpb+1
2440             ihpb(nhpb)=j
2441             jhpb(nhpb)=k
2442               dhpb(nhpb)=ddjk
2443             forcon(nhpb)=wfrag_(i)
2444             else if (constr_dist.eq.2) then
2445               if (ddjk.le.dist_cut) then
2446                 nhpb=nhpb+1
2447                 ihpb(nhpb)=j
2448                 jhpb(nhpb)=k
2449                 dhpb(nhpb)=ddjk
2450                 forcon(nhpb)=wfrag_(i)
2451               endif
2452             else
2453               nhpb=nhpb+1
2454               ihpb(nhpb)=j
2455               jhpb(nhpb)=k
2456               dhpb(nhpb)=ddjk
2457               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2458             endif
2459 #ifdef MPI
2460             if (.not.out1file .or. me.eq.king)
2461      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2462      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2463 #else
2464             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2465      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2466 #endif
2467           enddo
2468         enddo
2469         endif
2470       enddo
2471       do i=1,npair_
2472         if (wpair_(i).gt.0.0d0) then
2473         ii = ipair_(1,i)
2474         jj = ipair_(2,i)
2475         if (ii.gt.jj) then
2476           itemp=ii
2477           ii=jj
2478           jj=itemp
2479         endif
2480         do j=ifrag_(1,ii),ifrag_(2,ii)
2481           do k=ifrag_(1,jj),ifrag_(2,jj)
2482             nhpb=nhpb+1
2483             ihpb(nhpb)=j
2484             jhpb(nhpb)=k
2485             forcon(nhpb)=wpair_(i)
2486             dhpb(nhpb)=dist(j,k)
2487 #ifdef MPI
2488             if (.not.out1file .or. me.eq.king)
2489      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2490      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2491 #else
2492             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2493      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2494 #endif
2495           enddo
2496         enddo
2497         endif
2498       enddo
2499       print *,ndist_
2500       do i=1,ndist_
2501         if (constr_dist.eq.11) then
2502         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2503      &     ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2504         fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2505         else
2506 C        print *,"in else"
2507         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2508      &     ibecarb(i),forcon(nhpb+1)
2509         endif
2510         if (forcon(nhpb+1).gt.0.0d0) then
2511           nhpb=nhpb+1
2512           if (ibecarb(i).gt.0) then
2513             ihpb(i)=ihpb(i)+nres
2514             jhpb(i)=jhpb(i)+nres
2515           endif
2516           if (dhpb(nhpb).eq.0.0d0)
2517      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2518         endif
2519 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2520 C        if (forcon(nhpb+1).gt.0.0d0) then
2521 C          nhpb=nhpb+1
2522 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2523 #ifdef MPI
2524           if (.not.out1file .or. me.eq.king)
2525      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2526      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2527 #else
2528           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2529      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2530 #endif
2531
2532       enddo
2533       call flush(iout)
2534       return
2535       end
2536 c-------------------------------------------------------------------------------
2537 #ifdef WINIFL
2538       subroutine flush(iu)
2539       return
2540       end
2541 #endif
2542 #ifdef AIX
2543       subroutine flush(iu)
2544       call flush_(iu)
2545       return
2546       end
2547 #endif
2548 c------------------------------------------------------------------------------
2549       subroutine copy_to_tmp(source)
2550       include "DIMENSIONS"
2551       include "COMMON.IOUNITS"
2552       character*(*) source
2553       character* 256 tmpfile
2554       integer ilen
2555       external ilen
2556       logical ex
2557       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2558       inquire(file=tmpfile,exist=ex)
2559       if (ex) then
2560         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2561      &   " to temporary directory..."
2562         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2563         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2564       endif
2565       return
2566       end
2567 c------------------------------------------------------------------------------
2568       subroutine move_from_tmp(source)
2569       include "DIMENSIONS"
2570       include "COMMON.IOUNITS"
2571       character*(*) source
2572       integer ilen
2573       external ilen
2574       write (*,*) "Moving ",source(:ilen(source)),
2575      & " from temporary directory to working directory"
2576       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2577       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2578       return
2579       end
2580 c------------------------------------------------------------------------------
2581       subroutine random_init(seed)
2582 C
2583 C Initialize random number generator
2584 C
2585       implicit real*8 (a-h,o-z)
2586       include 'DIMENSIONS'
2587 #ifdef MPI
2588       include 'mpif.h'
2589       logical OKRandom, prng_restart
2590       real*8  r1
2591       integer iseed_array(4)
2592 #endif
2593       include 'COMMON.IOUNITS'
2594       include 'COMMON.TIME1'
2595       include 'COMMON.THREAD'
2596       include 'COMMON.SBRIDGE'
2597       include 'COMMON.CONTROL'
2598       include 'COMMON.MCM'
2599       include 'COMMON.MAP'
2600       include 'COMMON.HEADER'
2601       include 'COMMON.CSA'
2602       include 'COMMON.CHAIN'
2603       include 'COMMON.MUCA'
2604       include 'COMMON.MD'
2605       include 'COMMON.FFIELD'
2606       include 'COMMON.SETUP'
2607       iseed=-dint(dabs(seed))
2608       if (iseed.eq.0) then
2609         write (iout,'(/80(1h*)/20x,a/80(1h*))')
2610      &    'Random seed undefined. The program will stop.'
2611         write (*,'(/80(1h*)/20x,a/80(1h*))')
2612      &    'Random seed undefined. The program will stop.'
2613 #ifdef MPI
2614         call mpi_finalize(mpi_comm_world,ierr)
2615 #endif
2616         stop 'Bad random seed.'
2617       endif
2618 #ifdef MPI
2619       if (fg_rank.eq.0) then
2620       seed=seed*(me+1)+1
2621 #ifdef AMD64
2622       if(me.eq.king)
2623      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2624       OKRandom = prng_restart(me,iseed)
2625 #else
2626       do i=1,4
2627        tmp=65536.0d0**(4-i)
2628        iseed_array(i) = dint(seed/tmp)
2629        seed=seed-iseed_array(i)*tmp
2630       enddo
2631       if(me.eq.king)
2632      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2633      &                 (iseed_array(i),i=1,4)
2634       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2635      &                 (iseed_array(i),i=1,4)
2636       OKRandom = prng_restart(me,iseed_array)
2637 #endif
2638       if (OKRandom) then
2639 c        r1 = prng_next(me)
2640         r1=ran_number(0.0D0,1.0D0)
2641         if(me.eq.king)
2642      &   write (iout,*) 'ran_num',r1
2643         if (r1.lt.0.0d0) OKRandom=.false.
2644       endif
2645       if (.not.OKRandom) then
2646         write (iout,*) 'PRNG IS NOT WORKING!!!'
2647         print *,'PRNG IS NOT WORKING!!!'
2648         if (me.eq.0) then
2649          call flush(iout)
2650          call mpi_abort(mpi_comm_world,error_msg,ierr)
2651          stop
2652         else
2653          write (iout,*) 'too many processors for parallel prng'
2654          write (*,*) 'too many processors for parallel prng'
2655          call flush(iout)
2656          stop
2657         endif
2658       endif
2659       endif
2660 #else
2661       call vrndst(iseed)
2662       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
2663 #endif
2664       return
2665       end