source/unres/src-HCD-5D/readrtns_CSA.F

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