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