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