source/unres/src_MD-M/readrtns_CSA.F

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