2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.CHAIN'
13 C Read force-field parameters except weights
15 C Read job setup parameters
17 C Read control parameters for energy minimzation if required
18 if (minim) call read_minim
19 C Read MCM control parameters if required
20 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
21 C Read MD control parameters if reqjuired
22 if (modecalc.eq.12) call read_MDpar
23 C Read MREMD control parameters if required
24 if (modecalc.eq.14) then
28 C Read MUCA control parameters if required
29 if (lmuca) call read_muca
30 C Read CSA control parameters if required (from fort.40 if exists
31 C otherwise from general input file)
32 csa if (modecalc.eq.8) then
33 csa inquire (file="fort.40",exist=file_exist)
34 csa if (.not.file_exist) call csaread
36 cfmc if (modecalc.eq.10) call mcmfread
37 C Read molecule information, molecule geometry, energy-term weights, and
38 C restraints if requested
40 C Print restraint information
42 if (.not. out1file .or. me.eq.king) then
45 write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
46 & "The following",nhpb-nss,
47 & " distance restraints have been imposed:",
48 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
51 write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
52 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
59 c print *,"Processor",myrank," leaves READRTNS"
62 C-------------------------------------------------------------------------------
63 subroutine read_control
67 implicit real*8 (a-h,o-z)
71 logical OKRandom, prng_restart
74 include 'COMMON.IOUNITS'
75 include 'COMMON.TIME1'
76 include 'COMMON.THREAD'
77 include 'COMMON.SBRIDGE'
78 include 'COMMON.CONTROL'
81 include 'COMMON.HEADER'
82 csa include 'COMMON.CSA'
83 include 'COMMON.CHAIN'
86 include 'COMMON.FFIELD'
87 include 'COMMON.SETUP'
88 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
89 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
91 character*320 controlcard
97 read (INP,'(a)') titel
98 call card_concat(controlcard)
99 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
100 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
101 call reada(controlcard,'SEED',seed,0.0D0)
102 call random_init(seed)
103 C Set up the time limit (caution! The time must be input in minutes!)
104 read_cart=index(controlcard,'READ_CART').gt.0
105 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
106 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
107 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
108 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
109 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
110 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
111 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
112 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
113 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
114 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
115 call reada(controlcard,'DRMS',drms,0.1D0)
116 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
117 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
118 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
119 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
120 write (iout,'(a,f10.1)')'DRMS = ',drms
121 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
122 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
124 call readi(controlcard,'NZ_START',nz_start,0)
125 call readi(controlcard,'NZ_END',nz_end,0)
126 call readi(controlcard,'IZ_SC',iz_sc,0)
128 safety = 60.0d0*safety
131 call reada(controlcard,"T_BATH",t_bath,300.0d0)
132 minim=(index(controlcard,'MINIMIZE').gt.0)
133 dccart=(index(controlcard,'CART').gt.0)
134 overlapsc=(index(controlcard,'OVERLAP').gt.0)
135 overlapsc=.not.overlapsc
136 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
137 searchsc=.not.searchsc
138 sideadd=(index(controlcard,'SIDEADD').gt.0)
139 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
140 outpdb=(index(controlcard,'PDBOUT').gt.0)
141 outx=(index(controlcard,'XOUT').gt.0)
142 outmol2=(index(controlcard,'MOL2OUT').gt.0)
143 pdbref=(index(controlcard,'PDBREF').gt.0)
144 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
145 indpdb=index(controlcard,'PDBSTART')
146 extconf=(index(controlcard,'EXTCONF').gt.0)
147 call readi(controlcard,'IPRINT',iprint,0)
148 call readi(controlcard,'MAXGEN',maxgen,10000)
149 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
150 call readi(controlcard,"KDIAG",kdiag,0)
151 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
152 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
153 & write (iout,*) "RESCALE_MODE",rescale_mode
154 split_ene=index(controlcard,'SPLIT_ENE').gt.0
155 if (index(controlcard,'REGULAR').gt.0.0D0) then
156 call reada(controlcard,'WEIDIS',weidis,0.1D0)
160 call reada(controlcard,"CHECKGRAD_INC",checkgrad_inc,1.0d-4)
161 if (index(controlcard,'CHECKGRAD').gt.0) then
163 if (index(controlcard,'CART').gt.0) then
165 elseif (index(controlcard,'CARINT').gt.0) then
170 elseif (index(controlcard,'THREAD').gt.0) then
172 call readi(controlcard,'THREAD',nthread,0)
173 if (nthread.gt.0) then
174 call reada(controlcard,'WEIDIS',weidis,0.1D0)
177 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
178 stop 'Error termination in Read_Control.'
180 else if (index(controlcard,'MCMA').gt.0) then
182 else if (index(controlcard,'MCEE').gt.0) then
184 else if (index(controlcard,'MULTCONF').gt.0) then
186 else if (index(controlcard,'MAP').gt.0) then
188 call readi(controlcard,'MAP',nmap,0)
189 else if (index(controlcard,'CSA').gt.0) then
190 write(*,*) "CSA not supported in this version"
193 crc else if (index(controlcard,'ZSCORE').gt.0) then
195 crc ZSCORE is rm from UNRES, modecalc=9 is available
198 cfcm else if (index(controlcard,'MCMF').gt.0) then
200 else if (index(controlcard,'SOFTREG').gt.0) then
202 else if (index(controlcard,'CHECK_BOND').gt.0) then
204 else if (index(controlcard,'TEST').gt.0) then
206 else if (index(controlcard,'MD').gt.0) then
208 else if (index(controlcard,'RE ').gt.0) then
212 lmuca=index(controlcard,'MUCA').gt.0
213 call readi(controlcard,'MUCADYN',mucadyn,0)
214 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
215 if (lmuca .and. (me.eq.king .or. .not.out1file ))
217 write (iout,*) 'MUCADYN=',mucadyn
218 write (iout,*) 'MUCASMOOTH=',muca_smooth
221 iscode=index(controlcard,'ONE_LETTER')
222 indphi=index(controlcard,'PHI')
223 indback=index(controlcard,'BACK')
224 iranconf=index(controlcard,'RAND_CONF')
225 i2ndstr=index(controlcard,'USE_SEC_PRED')
226 gradout=index(controlcard,'GRADOUT').gt.0
227 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
229 if(me.eq.king.or..not.out1file)
230 & write (iout,'(2a)') diagmeth(kdiag),
231 & ' routine used to diagonalize matrices.'
234 c--------------------------------------------------------------------------
235 subroutine read_REMDpar
239 implicit real*8 (a-h,o-z)
241 include 'COMMON.IOUNITS'
242 include 'COMMON.TIME1'
245 include 'COMMON.LANGEVIN'
247 include 'COMMON.LANGEVIN.lang0'
249 include 'COMMON.INTERACT'
250 include 'COMMON.NAMES'
252 include 'COMMON.REMD'
253 include 'COMMON.CONTROL'
254 include 'COMMON.SETUP'
256 character*320 controlcard
257 character*3200 controlcard1
258 integer iremd_m_total
260 if(me.eq.king.or..not.out1file)
261 & write (iout,*) "REMD setup"
263 call card_concat(controlcard)
264 call readi(controlcard,"NREP",nrep,3)
265 call readi(controlcard,"NSTEX",nstex,1000)
266 call reada(controlcard,"RETMIN",retmin,10.0d0)
267 call reada(controlcard,"RETMAX",retmax,1000.0d0)
268 mremdsync=(index(controlcard,'SYNC').gt.0)
269 call readi(controlcard,"NSYN",i_sync_step,100)
270 restart1file=(index(controlcard,'REST1FILE').gt.0)
271 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
272 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
273 if(max_cache_traj_use.gt.max_cache_traj)
274 & max_cache_traj_use=max_cache_traj
275 if(me.eq.king.or..not.out1file) then
276 cd if (traj1file) then
277 crc caching is in testing - NTWX is not ignored
278 cd write (iout,*) "NTWX value is ignored"
279 cd write (iout,*) " trajectory is stored to one file by master"
280 cd write (iout,*) " before exchange at NSTEX intervals"
282 write (iout,*) "NREP= ",nrep
283 write (iout,*) "NSTEX= ",nstex
284 write (iout,*) "SYNC= ",mremdsync
285 write (iout,*) "NSYN= ",i_sync_step
286 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
289 t_exchange_only=(index(controlcard,'TONLY').gt.0)
290 call readi(controlcard,"HREMD",hremd,0)
291 if((me.eq.king.or..not.out1file).and.hremd.gt.0) then
292 write (iout,*) "Hamiltonian REMD with ",hremd," sets of weights"
294 if(usampl.and.hremd.gt.0) then
296 & "========== ERROR: USAMPL and HREMD cannot be used together"
298 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
305 if (index(controlcard,'TLIST').gt.0) then
307 call card_concat(controlcard1)
308 read(controlcard1,*) (remd_t(i),i=1,nrep)
309 if(me.eq.king.or..not.out1file)
310 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
313 if (index(controlcard,'MLIST').gt.0) then
315 call card_concat(controlcard1)
316 read(controlcard1,*) (remd_m(i),i=1,nrep)
317 if(me.eq.king.or..not.out1file) then
318 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
321 iremd_m_total=iremd_m_total+remd_m(i)
324 write (iout,*) 'Total number of replicas ',
325 & iremd_m_total*hremd
327 write (iout,*) 'Total number of replicas ',iremd_m_total
331 if(me.eq.king.or..not.out1file)
332 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
335 c--------------------------------------------------------------------------
336 subroutine read_MDpar
340 implicit real*8 (a-h,o-z)
342 include 'COMMON.IOUNITS'
343 include 'COMMON.TIME1'
346 include 'COMMON.LANGEVIN'
348 include 'COMMON.LANGEVIN.lang0'
350 include 'COMMON.INTERACT'
351 include 'COMMON.NAMES'
353 include 'COMMON.SETUP'
354 include 'COMMON.CONTROL'
355 include 'COMMON.SPLITELE'
357 character*320 controlcard
359 call card_concat(controlcard)
360 call readi(controlcard,"NSTEP",n_timestep,1000000)
361 call readi(controlcard,"NTWE",ntwe,100)
362 call readi(controlcard,"NTWX",ntwx,1000)
363 call reada(controlcard,"DT",d_time,1.0d-1)
364 call reada(controlcard,"DVMAX",dvmax,2.0d1)
365 call reada(controlcard,"DAMAX",damax,1.0d1)
366 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
367 call readi(controlcard,"LANG",lang,0)
368 RESPA = index(controlcard,"RESPA") .gt. 0
369 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
370 ntime_split0=ntime_split
371 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
372 ntime_split0=ntime_split
373 call reada(controlcard,"R_CUT",r_cut,2.0d0)
374 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
375 rest = index(controlcard,"REST").gt.0
376 tbf = index(controlcard,"TBF").gt.0
377 call readi(controlcard,"HMC",hmc,0)
378 tnp = index(controlcard,"NOSEPOINCARE99").gt.0
379 tnp1 = index(controlcard,"NOSEPOINCARE01").gt.0
380 tnh = index(controlcard,"NOSEHOOVER96").gt.0
381 if (RESPA.and.tnh)then
382 xiresp = index(controlcard,"XIRESP").gt.0
384 call reada(controlcard,"Q_NP",Q_np,0.1d0)
385 usampl = index(controlcard,"USAMPL").gt.0
386 mdpdb = index(controlcard,"MDPDB").gt.0
387 call reada(controlcard,"T_BATH",t_bath,300.0d0)
388 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
389 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
390 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
391 if (count_reset_moment.eq.0) count_reset_moment=1000000000
392 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
393 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
394 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
395 if (count_reset_vel.eq.0) count_reset_vel=1000000000
396 large = index(controlcard,"LARGE").gt.0
397 print_compon = index(controlcard,"PRINT_COMPON").gt.0
398 rattle = index(controlcard,"RATTLE").gt.0
399 preminim = index(controlcard,"PREMINIM").gt.0
401 dccart=(index(controlcard,'CART').gt.0)
404 c if performing umbrella sampling, fragments constrained are read from the fragment file
410 if(me.eq.king.or..not.out1file) then
412 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
414 write (iout,'(a)') "The units are:"
415 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
416 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
417 & " acceleration: angstrom/(48.9 fs)**2"
418 write (iout,'(a)') "energy: kcal/mol, temperature: K"
420 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
421 write (iout,'(a60,f10.5,a)')
422 & "Initial time step of numerical integration:",d_time,
424 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
426 write (iout,'(2a,i4,a)')
427 & "A-MTS algorithm used; initial time step for fast-varying",
428 & " short-range forces split into",ntime_split," steps."
429 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
430 & r_cut," lambda",rlamb
432 write (iout,'(2a,f10.5)')
433 & "Maximum acceleration threshold to reduce the time step",
434 & "/increase split number:",damax
435 write (iout,'(2a,f10.5)')
436 & "Maximum predicted energy drift to reduce the timestep",
437 & "/increase split number:",edriftmax
438 write (iout,'(a60,f10.5)')
439 & "Maximum velocity threshold to reduce velocities:",dvmax
440 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
441 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
442 if (rattle) write (iout,'(a60)')
443 & "Rattle algorithm used to constrain the virtual bonds"
444 if (preminim .or. iranconf.gt.0) then
446 & "Initial structure will be energy-minimized"
451 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
452 call reada(controlcard,"RWAT",rwat,1.4d0)
453 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
454 surfarea=index(controlcard,"SURFAREA").gt.0
455 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
456 if(me.eq.king.or..not.out1file)then
457 write (iout,'(/a,$)') "Langevin dynamics calculation"
460 & " with direct integration of Langevin equations"
461 else if (lang.eq.2) then
462 write (iout,'(a/)') " with TINKER stochasic MD integrator"
463 else if (lang.eq.3) then
464 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
465 else if (lang.eq.4) then
466 write (iout,'(a/)') " in overdamped mode"
468 write (iout,'(//a,i5)')
469 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
472 write (iout,'(a60,f10.5)') "Temperature:",t_bath
473 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
474 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
475 write (iout,'(a60,f10.5)')
476 & "Scaling factor of the friction forces:",scal_fric
477 if (surfarea) write (iout,'(2a,i10,a)')
478 & "Friction coefficients will be scaled by solvent-accessible",
479 & " surface area every",reset_fricmat," steps."
481 c Calculate friction coefficients and bounds of stochastic forces
482 eta=6*pi*cPoise*etawat
483 if(me.eq.king.or..not.out1file)
484 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
486 gamp=scal_fric*(pstok+rwat)*eta
487 stdfp=dsqrt(2*Rb*t_bath/d_time)
489 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
490 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
492 if(me.eq.king.or..not.out1file)then
493 write (iout,'(/2a/)')
494 & "Radii of site types and friction coefficients and std's of",
495 & " stochastic forces of fully exposed sites"
496 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
498 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
499 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
503 if(me.eq.king.or..not.out1file)then
504 write (iout,'(a)') "Berendsen bath calculation"
505 write (iout,'(a60,f10.5)') "Temperature:",t_bath
506 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
508 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
509 & count_reset_moment," steps"
511 & write (iout,'(a,i10,a)')
512 & "Velocities will be reset at random every",count_reset_vel,
515 else if (tnp .or. tnp1 .or. tnh) then
516 if (tnp .or. tnp1) then
517 write (iout,'(a)') "Nose-Poincare bath calculation"
518 if (tnp) write (iout,'(a)')
519 & "J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird"
520 if (tnp1) write (iout,'(a)') "JPSJ 70 75 (2001) S. Nose"
522 write (iout,'(a)') "Nose-Hoover bath calculation"
523 write (iout,'(a)') "Mol.Phys. 87 1117 (1996) Martyna et al."
533 WDTI(i) = 1.0*d_time/nresn
540 write (iout,'(a)') "NVT-XI-RESPA algorithm"
542 write (iout,'(a)') "NVT-XO-RESPA algorithm"
545 WDTIi(i) = 1.0*d_time/nresn/ntime_split
553 write (iout,'(a60,f10.5)') "Temperature:",t_bath
554 write (iout,'(a60,f10.5)') "Q =",Q_np
556 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
557 & count_reset_moment," steps"
559 & write (iout,'(a,i10,a)')
560 & "Velocities will be reset at random every",count_reset_vel,
563 else if (hmc.gt.0) then
564 write (iout,'(a)') "Hybrid Monte Carlo calculation"
565 write (iout,'(a60,f10.5)') "Temperature:",t_bath
566 write (iout,'(a60,i10)')
567 & "Number of MD steps between Metropolis tests:",hmc
570 if(me.eq.king.or..not.out1file)
571 & write (iout,'(a31)') "Microcanonical mode calculation"
573 if(me.eq.king.or..not.out1file)then
574 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
576 write(iout,*) "MD running with constraints."
577 write(iout,*) "Equilibration time ", eq_time, " mtus."
578 write(iout,*) "Constraining ", nfrag," fragments."
579 write(iout,*) "Length of each fragment, weight and q0:"
581 write (iout,*) "Set of restraints #",iset
583 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
584 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
586 write(iout,*) "constraints between ", npair, "fragments."
587 write(iout,*) "constraint pairs, weights and q0:"
589 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
590 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
592 write(iout,*) "angle constraints within ", nfrag_back,
593 & "backbone fragments."
594 write(iout,*) "fragment, weights:"
596 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
597 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
598 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
601 iset=mod(kolor,nset)+1
604 if(me.eq.king.or..not.out1file)
605 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
608 c------------------------------------------------------------------------------
611 C Read molecular data.
613 implicit real*8 (a-h,o-z)
619 include 'COMMON.IOUNITS'
622 include 'COMMON.INTERACT'
623 include 'COMMON.LOCAL'
624 include 'COMMON.NAMES'
625 include 'COMMON.CHAIN'
626 include 'COMMON.FFIELD'
627 include 'COMMON.SBRIDGE'
628 include 'COMMON.HEADER'
629 include 'COMMON.CONTROL'
630 include 'COMMON.DBASE'
631 include 'COMMON.THREAD'
632 include 'COMMON.CONTACTS'
633 include 'COMMON.TORCNSTR'
634 include 'COMMON.TIME1'
635 include 'COMMON.BOUNDS'
637 include 'COMMON.REMD'
638 include 'COMMON.SETUP'
639 character*4 sequence(maxres)
641 double precision x(maxvar)
642 character*256 pdbfile
643 character*320 weightcard
644 character*80 weightcard_t,ucase
645 dimension itype_pdb(maxres)
646 common /pizda/ itype_pdb
647 logical seq_comp,fail
648 double precision energia(0:n_ene)
655 C Read weights of the subsequent energy terms.
668 if(me.eq.king.or..not.out1file) then
669 write (iout,*) 'Reading ',hremd,' sets of weights for HREMD'
670 write (iout,*) 'Current weights for processor ',
671 & me,' set ',i2set(me)
675 call card_concat(weightcard)
676 call reada(weightcard,'WLONG',wlong,1.0D0)
677 call reada(weightcard,'WSC',wsc,wlong)
678 call reada(weightcard,'WSCP',wscp,wlong)
679 call reada(weightcard,'WELEC',welec,1.0D0)
680 call reada(weightcard,'WVDWPP',wvdwpp,welec)
681 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
682 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
683 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
684 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
685 call reada(weightcard,'WTURN3',wturn3,1.0D0)
686 call reada(weightcard,'WTURN4',wturn4,1.0D0)
687 call reada(weightcard,'WTURN6',wturn6,1.0D0)
688 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
689 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
690 call reada(weightcard,'WBOND',wbond,1.0D0)
691 call reada(weightcard,'WTOR',wtor,1.0D0)
692 call reada(weightcard,'WTORD',wtor_d,1.0D0)
693 call reada(weightcard,'WANG',wang,1.0D0)
694 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
695 call reada(weightcard,'SCAL14',scal14,0.4D0)
696 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
697 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
698 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
699 call reada(weightcard,'TEMP0',temp0,300.0d0)
700 if (index(weightcard,'SOFT').gt.0) ipot=6
701 C 12/1/95 Added weight for the multi-body term WCORR
702 call reada(weightcard,'WCORRH',wcorr,1.0D0)
703 if (wcorr4.gt.0.0d0) wcorr=wcorr4
711 hweights(i,7)=wel_loc
714 hweights(i,10)=wturn6
716 hweights(i,12)=wscloc
718 hweights(i,14)=wtor_d
719 hweights(i,15)=wstrain
720 hweights(i,16)=wvdwpp
722 hweights(i,18)=scal14
723 hweights(i,21)=wsccor
728 weights(i)=hweights(i2set(me),i)
752 call card_concat(weightcard)
753 call reada(weightcard,'WLONG',wlong,1.0D0)
754 call reada(weightcard,'WSC',wsc,wlong)
755 call reada(weightcard,'WSCP',wscp,wlong)
756 call reada(weightcard,'WELEC',welec,1.0D0)
757 call reada(weightcard,'WVDWPP',wvdwpp,welec)
758 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
759 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
760 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
761 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
762 call reada(weightcard,'WTURN3',wturn3,1.0D0)
763 call reada(weightcard,'WTURN4',wturn4,1.0D0)
764 call reada(weightcard,'WTURN6',wturn6,1.0D0)
765 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
766 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
767 call reada(weightcard,'WBOND',wbond,1.0D0)
768 call reada(weightcard,'WTOR',wtor,1.0D0)
769 call reada(weightcard,'WTORD',wtor_d,1.0D0)
770 call reada(weightcard,'WANG',wang,1.0D0)
771 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
772 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
773 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
774 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
775 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
776 call reada(weightcard,'SCAL14',scal14,0.4D0)
777 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
778 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
779 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
780 call reada(weightcard,'TEMP0',temp0,300.0d0)
781 if (index(weightcard,'SOFT').gt.0) ipot=6
782 C 12/1/95 Added weight for the multi-body term WCORR
783 call reada(weightcard,'WCORRH',wcorr,1.0D0)
784 if (wcorr4.gt.0.0d0) wcorr=wcorr4
805 weights(25)=wdfa_dist
808 weights(28)=wdfa_beta
810 if(me.eq.king.or..not.out1file)
811 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
812 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
814 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
815 10 format (/'Energy-term weights (unscaled):'//
816 & 'WSCC= ',f10.6,' (SC-SC)'/
817 & 'WSCP= ',f10.6,' (SC-p)'/
818 & 'WELEC= ',f10.6,' (p-p electr)'/
819 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
820 & 'WBOND= ',f10.6,' (stretching)'/
821 & 'WANG= ',f10.6,' (bending)'/
822 & 'WSCLOC= ',f10.6,' (SC local)'/
823 & 'WTOR= ',f10.6,' (torsional)'/
824 & 'WTORD= ',f10.6,' (double torsional)'/
825 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
826 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
827 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
828 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
829 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
830 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
831 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
832 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
833 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
834 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
835 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
836 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
837 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
838 if(me.eq.king.or..not.out1file)then
839 if (wcorr4.gt.0.0d0) then
840 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
841 & 'between contact pairs of peptide groups'
842 write (iout,'(2(a,f5.3/))')
843 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
844 & 'Range of quenching the correlation terms:',2*delt_corr
845 else if (wcorr.gt.0.0d0) then
846 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
847 & 'between contact pairs of peptide groups'
849 write (iout,'(a,f8.3)')
850 & 'Scaling factor of 1,4 SC-p interactions:',scal14
851 write (iout,'(a,f8.3)')
852 & 'General scaling factor of SC-p interactions:',scalscp
854 r0_corr=cutoff_corr-delt_corr
856 aad(i,1)=scalscp*aad(i,1)
857 aad(i,2)=scalscp*aad(i,2)
858 bad(i,1)=scalscp*bad(i,1)
859 bad(i,2)=scalscp*bad(i,2)
861 call rescale_weights(t_bath)
862 if(me.eq.king.or..not.out1file)
863 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
864 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
866 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
867 22 format (/'Energy-term weights (scaled):'//
868 & 'WSCC= ',f10.6,' (SC-SC)'/
869 & 'WSCP= ',f10.6,' (SC-p)'/
870 & 'WELEC= ',f10.6,' (p-p electr)'/
871 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
872 & 'WBOND= ',f10.6,' (stretching)'/
873 & 'WANG= ',f10.6,' (bending)'/
874 & 'WSCLOC= ',f10.6,' (SC local)'/
875 & 'WTOR= ',f10.6,' (torsional)'/
876 & 'WTORD= ',f10.6,' (double torsional)'/
877 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
878 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
879 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
880 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
881 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
882 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
883 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
884 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
885 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
886 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
887 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
888 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
889 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
890 if(me.eq.king.or..not.out1file)
891 & write (iout,*) "Reference temperature for weights calculation:",
893 call reada(weightcard,"D0CM",d0cm,3.78d0)
894 call reada(weightcard,"AKCM",akcm,15.1d0)
895 call reada(weightcard,"AKTH",akth,11.0d0)
896 call reada(weightcard,"AKCT",akct,12.0d0)
897 call reada(weightcard,"V1SS",v1ss,-1.08d0)
898 call reada(weightcard,"V2SS",v2ss,7.61d0)
899 call reada(weightcard,"V3SS",v3ss,13.7d0)
900 call reada(weightcard,"EBR",ebr,-5.50D0)
901 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
903 dyn_ss_mask(i)=.false.
907 dyn_ssbond_ij(i,j)=1.0d300
910 call reada(weightcard,"HT",Ht,0.0D0)
912 ss_depth=ebr/wsc-0.25*eps(1,1)
913 Ht=Ht/wsc-0.25*eps(1,1)
914 akcm=akcm*wstrain/wsc
915 akth=akth*wstrain/wsc
916 akct=akct*wstrain/wsc
917 v1ss=v1ss*wstrain/wsc
918 v2ss=v2ss*wstrain/wsc
919 v3ss=v3ss*wstrain/wsc
921 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
924 if(me.eq.king.or..not.out1file) then
925 write (iout,*) "Parameters of the SS-bond potential:"
926 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
928 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
929 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
930 write (iout,*)" HT",Ht
931 print *,'indpdb=',indpdb,' pdbref=',pdbref
933 if (indpdb.gt.0 .or. pdbref) then
934 read(inp,'(a)') pdbfile
935 if(me.eq.king.or..not.out1file)
936 & write (iout,'(2a)') 'PDB data will be read from file ',
937 & pdbfile(:ilen(pdbfile))
938 open(ipdbin,file=pdbfile,status='old',err=33)
940 33 write (iout,'(a)') 'Error opening PDB file.'
943 c print *,'Begin reading pdb data'
952 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
953 & (crefjlee(j,i+nres),j=1,3)
956 c print *,'Finished reading pdb data'
957 if(me.eq.king.or..not.out1file)
958 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
959 & ' nstart_sup=',nstart_sup
961 itype_pdb(i)=itype(i)
965 nct=nstart_sup+nsup-1
966 call contact(.false.,ncont_ref,icont_ref,co)
969 C Following 2 lines for diagnostics; comment out if not needed
970 write (iout,*) "Before sideadd"
972 if(me.eq.king.or..not.out1file)
973 & write(iout,*)'Adding sidechains'
980 do while (fail.and.nsi.le.maxsi)
981 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
984 if(fail) write(iout,*)'Adding sidechain failed for res ',
985 & i,' after ',nsi,' trials'
988 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
991 C Following 2 lines for diagnostics; comment out if not needed
992 c write (iout,*) "After sideadd"
995 if (indpdb.eq.0) then
996 C Read sequence if not taken from the pdb file.
998 c print *,'nres=',nres
999 if (iscode.gt.0) then
1000 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
1002 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
1004 C Convert sequence to numeric code
1006 itype(i)=rescode(i,sequence(i),iscode)
1008 if (itype(2).eq.10.and.itype(1).eq.ntyp1) then
1010 & "Glycine is the first full residue, initial dummy deleted"
1016 if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then
1018 & "Glycine is the last full residue, terminal dummy deleted"
1022 C Assign initial virtual bond lengths
1028 vbld(i+nres)=dsc(itype(i))
1029 vbld_inv(i+nres)=dsc_inv(itype(i))
1030 c write (iout,*) "i",i," itype",itype(i),
1031 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
1035 c print '(20i4)',(itype(i),i=1,nres)
1038 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
1040 if (itype(i).eq.21) then
1044 else if (itype(i+1).ne.20) then
1046 else if (itype(i).ne.20) then
1053 if(me.eq.king.or..not.out1file)then
1054 write (iout,*) "ITEL"
1056 write (iout,*) i,itype(i),itel(i)
1058 print *,'Call Read_Bridge.'
1061 C 8/13/98 Set limits to generating the dihedral angles
1066 read (inp,*) ndih_constr
1067 if (ndih_constr.gt.0) then
1069 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1070 if(me.eq.king.or..not.out1file)then
1072 & 'There are',ndih_constr,' constraints on phi angles.'
1074 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1078 phi0(i)=deg2rad*phi0(i)
1079 drange(i)=deg2rad*drange(i)
1081 if(me.eq.king.or..not.out1file)
1082 & write (iout,*) 'FTORS',ftors
1085 phibound(1,ii) = phi0(i)-drange(i)
1086 phibound(2,ii) = phi0(i)+drange(i)
1091 if (me.eq.king) then
1093 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1095 write (iout,'(a3,i5,2f10.1)')
1096 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1102 cd print *,'NNT=',NNT,' NCT=',NCT
1103 if (itype(1).eq.21) nnt=2
1104 if (itype(nres).eq.21) nct=nct-1
1106 C Bartek:READ init_vars
1107 C Initialize variables!
1108 C Juyong:READ read_info
1109 C READ fragment information!!
1110 C both routines should be in dfa.F file!!
1113 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
1114 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
1116 print*, 'init_dfa_vars finished!'
1118 print*, 'read_dfa_info finished!'
1123 if(me.eq.king.or..not.out1file)
1124 & write (iout,'(a,i3)') 'nsup=',nsup
1126 if (nsup.le.(nct-nnt+1)) then
1127 do i=0,nct-nnt+1-nsup
1128 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1134 & 'Error - sequences to be superposed do not match.'
1137 do i=0,nsup-(nct-nnt+1)
1138 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1140 nstart_sup=nstart_sup+i
1146 & 'Error - sequences to be superposed do not match.'
1149 if (nsup.eq.0) nsup=nct-nnt
1150 if (nstart_sup.eq.0) nstart_sup=nnt
1151 if (nstart_seq.eq.0) nstart_seq=nnt
1152 if(me.eq.king.or..not.out1file)
1153 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1154 & ' nstart_seq=',nstart_seq
1156 c--- Zscore rms -------
1157 if (nz_start.eq.0) nz_start=nnt
1158 if (nz_end.eq.0 .and. nsup.gt.0) then
1160 else if (nz_end.eq.0) then
1163 if(me.eq.king.or..not.out1file)then
1164 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1165 write (iout,*) 'IZ_SC=',iz_sc
1167 c----------------------
1170 if (.not.pdbref) then
1171 call read_angles(inp,*38)
1173 38 write (iout,'(a)') 'Error reading reference structure.'
1175 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1176 stop 'Error reading reference structure'
1180 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1189 call contact(.true.,ncont_ref,icont_ref,co)
1191 if(me.eq.king.or..not.out1file)
1192 & write (iout,*) 'Contact order:',co
1194 if(me.eq.king.or..not.out1file)
1195 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1198 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1200 if(me.eq.king.or..not.out1file)
1201 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1202 & icont_ref(1,i),' ',
1203 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1207 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1208 if (constr_dist.gt.0) then
1209 call read_dist_constr
1213 if (constr_homology.gt.0) then
1214 call read_constr_homology
1215 if (indpdb.gt.0 .or. pdbref) then
1218 c(j,i)=crefjlee(j,i)
1219 cref(j,i)=crefjlee(j,i)
1224 write (iout,*) "Array C"
1226 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1227 & (c(j,i+nres),j=1,3)
1229 write (iout,*) "Array Cref"
1231 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1232 & (cref(j,i+nres),j=1,3)
1235 call int_from_cart1(.false.)
1236 call sc_loc_geom(.false.)
1238 thetaref(i)=theta(i)
1243 dc(j,i)=c(j,i+1)-c(j,i)
1244 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1249 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1250 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1258 if (nhpb.gt.0) call hpb_partition
1259 c write (iout,*) "After read_dist_constr nhpb",nhpb
1261 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1262 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1263 & modecalc.ne.10) then
1264 C If input structure hasn't been supplied from the PDB file read or generate
1266 if (iranconf.eq.0 .and. .not. extconf) then
1267 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1268 & write (iout,'(a)') 'Initial geometry will be read in.'
1270 c read (inp,*) time,potE,uconst,t_bath,
1271 c & nss,(ihpb(j),jhpb(j),j=1,nss), nn, (qfrag(i),i=1,nn)
1272 read(inp,'(8f10.5)',end=36,err=36)
1273 & ((c(l,k),l=1,3),k=1,nres),
1274 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1275 call int_from_cart1(.false.)
1278 dc(j,i)=c(j,i+1)-c(j,i)
1279 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1283 if (itype(i).ne.10) then
1285 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1286 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1292 call read_angles(inp,*36)
1295 36 write (iout,'(a)') 'Error reading angle file.'
1297 call mpi_finalize( MPI_COMM_WORLD,IERR )
1299 stop 'Error reading angle file.'
1301 else if (extconf) then
1302 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1303 & write (iout,'(a)') 'Extended chain initial geometry.'
1305 theta(i)=90d0*deg2rad
1308 phi(i)=180d0*deg2rad
1311 alph(i)=110d0*deg2rad
1314 omeg(i)=-120d0*deg2rad
1317 if(me.eq.king.or..not.out1file)
1318 & write (iout,'(a)') 'Random-generated initial geometry.'
1322 if (me.eq.king .or. fg_rank.eq.0 .and. (
1323 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1327 call gen_rand_conf(itmp,*30)
1329 30 write (iout,*) 'Failed to generate random conformation',
1330 & ', itrial=',itrial
1331 write (*,*) 'Processor:',me,
1332 & ' Failed to generate random conformation',
1342 write (iout,'(a,i3,a)') 'Processor:',me,
1343 & ' error in generating random conformation.'
1344 write (*,'(a,i3,a)') 'Processor:',me,
1345 & ' error in generating random conformation.'
1348 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1355 elseif (modecalc.eq.4) then
1356 read (inp,'(a)') intinname
1357 open (intin,file=intinname,status='old',err=333)
1358 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1359 & write (iout,'(a)') 'intinname',intinname
1360 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1362 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1364 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1366 stop 'Error opening angle file.'
1370 C Generate distance constraints, if the PDB structure is to be regularized.
1371 if (nthread.gt.0) then
1372 call read_threadbase
1374 write (iout,*) "READRTNS: Calling setup_var"
1377 if (me.eq.king .or. .not. out1file)
1379 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1380 write (iout,'(/a,i3,a)')
1381 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1382 write (iout,'(20i4)') (iss(i),i=1,ns)
1384 write(iout,*)"Running with dynamic disulfide-bond formation"
1386 write (iout,'(/a/)') 'Pre-formed links are:'
1392 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1393 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1399 if (ns.gt.0.and.dyn_ss) then
1403 forcon(i-nss)=forcon(i)
1410 dyn_ss_mask(iss(i))=.true.
1413 if (i2ndstr.gt.0) call secstrp2dihc
1414 c call geom_to_var(nvar,x)
1415 c call etotal(energia(0))
1416 c call enerprint(energia(0))
1417 c call briefout(0,etot)
1419 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1420 cd write (iout,'(a)') 'Variable list:'
1421 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1423 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1424 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1425 & 'Processor',myrank,': end reading molecular data.'
1429 c--------------------------------------------------------------------------
1430 logical function seq_comp(itypea,itypeb,length)
1432 integer length,itypea(length),itypeb(length)
1435 if (itypea(i).ne.itypeb(i)) then
1443 c-----------------------------------------------------------------------------
1444 subroutine read_bridge
1445 C Read information about disulfide bridges.
1446 implicit real*8 (a-h,o-z)
1447 include 'DIMENSIONS'
1451 include 'COMMON.IOUNITS'
1452 include 'COMMON.GEO'
1453 include 'COMMON.VAR'
1454 include 'COMMON.INTERACT'
1455 include 'COMMON.LOCAL'
1456 include 'COMMON.NAMES'
1457 include 'COMMON.CHAIN'
1458 include 'COMMON.FFIELD'
1459 include 'COMMON.SBRIDGE'
1460 include 'COMMON.HEADER'
1461 include 'COMMON.CONTROL'
1462 include 'COMMON.DBASE'
1463 include 'COMMON.THREAD'
1464 include 'COMMON.TIME1'
1465 include 'COMMON.SETUP'
1466 C Read bridging residues.
1467 read (inp,*) ns,(iss(i),i=1,ns)
1469 if(me.eq.king.or..not.out1file)
1470 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1471 C Check whether the specified bridging residues are cystines.
1473 if (itype(iss(i)).ne.1) then
1474 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1475 & 'Do you REALLY think that the residue ',
1476 & restyp(itype(iss(i))),i,
1477 & ' can form a disulfide bridge?!!!'
1478 write (*,'(2a,i3,a)')
1479 & 'Do you REALLY think that the residue ',
1480 & restyp(itype(iss(i))),i,
1481 & ' can form a disulfide bridge?!!!'
1483 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1488 C Read preformed bridges.
1490 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1492 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1495 C Check if the residues involved in bridges are in the specified list of
1496 C bridging residues.
1499 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1500 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1501 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1502 & ' contains residues present in other pairs.'
1503 write (*,'(a,i3,a)') 'Disulfide pair',i,
1504 & ' contains residues present in other pairs.'
1506 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1512 if (ihpb(i).eq.iss(j)) goto 10
1514 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1517 if (jhpb(i).eq.iss(j)) goto 20
1519 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1525 ihpb(i)=ihpb(i)+nres
1526 jhpb(i)=jhpb(i)+nres
1532 c----------------------------------------------------------------------------
1533 subroutine read_x(kanal,*)
1534 implicit real*8 (a-h,o-z)
1535 include 'DIMENSIONS'
1536 include 'COMMON.GEO'
1537 include 'COMMON.VAR'
1538 include 'COMMON.CHAIN'
1539 include 'COMMON.IOUNITS'
1540 include 'COMMON.CONTROL'
1541 include 'COMMON.LOCAL'
1542 include 'COMMON.INTERACT'
1543 c Read coordinates from input
1545 read(kanal,'(8f10.5)',end=10,err=10)
1546 & ((c(l,k),l=1,3),k=1,nres),
1547 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1550 c(j,2*nres)=c(j,nres)
1552 call int_from_cart1(.false.)
1555 dc(j,i)=c(j,i+1)-c(j,i)
1556 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1560 if (itype(i).ne.10) then
1562 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1563 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1571 c----------------------------------------------------------------------------
1572 subroutine read_threadbase
1573 implicit real*8 (a-h,o-z)
1574 include 'DIMENSIONS'
1575 include 'COMMON.IOUNITS'
1576 include 'COMMON.GEO'
1577 include 'COMMON.VAR'
1578 include 'COMMON.INTERACT'
1579 include 'COMMON.LOCAL'
1580 include 'COMMON.NAMES'
1581 include 'COMMON.CHAIN'
1582 include 'COMMON.FFIELD'
1583 include 'COMMON.SBRIDGE'
1584 include 'COMMON.HEADER'
1585 include 'COMMON.CONTROL'
1586 include 'COMMON.DBASE'
1587 include 'COMMON.THREAD'
1588 include 'COMMON.TIME1'
1589 C Read pattern database for threading.
1590 read (icbase,*) nseq
1592 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1593 & nres_base(2,i),nres_base(3,i)
1594 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1596 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1597 c & nres_base(2,i),nres_base(3,i)
1598 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1602 if (weidis.eq.0.0D0) weidis=0.1D0
1611 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1612 write (iout,'(a,i5)') 'nexcl: ',nexcl
1613 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1616 c------------------------------------------------------------------------------
1617 subroutine setup_var
1618 implicit real*8 (a-h,o-z)
1619 include 'DIMENSIONS'
1620 include 'COMMON.IOUNITS'
1621 include 'COMMON.GEO'
1622 include 'COMMON.VAR'
1623 include 'COMMON.INTERACT'
1624 include 'COMMON.LOCAL'
1625 include 'COMMON.NAMES'
1626 include 'COMMON.CHAIN'
1627 include 'COMMON.FFIELD'
1628 include 'COMMON.SBRIDGE'
1629 include 'COMMON.HEADER'
1630 include 'COMMON.CONTROL'
1631 include 'COMMON.DBASE'
1632 include 'COMMON.THREAD'
1633 include 'COMMON.TIME1'
1634 C Set up variable list.
1640 if (itype(i).ne.10) then
1642 ialph(i,1)=nvar+nside
1646 if (indphi.gt.0) then
1648 else if (indback.gt.0) then
1653 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1656 c----------------------------------------------------------------------------
1657 subroutine gen_dist_constr
1658 C Generate CA distance constraints.
1659 implicit real*8 (a-h,o-z)
1660 include 'DIMENSIONS'
1661 include 'COMMON.IOUNITS'
1662 include 'COMMON.GEO'
1663 include 'COMMON.VAR'
1664 include 'COMMON.INTERACT'
1665 include 'COMMON.LOCAL'
1666 include 'COMMON.NAMES'
1667 include 'COMMON.CHAIN'
1668 include 'COMMON.FFIELD'
1669 include 'COMMON.SBRIDGE'
1670 include 'COMMON.HEADER'
1671 include 'COMMON.CONTROL'
1672 include 'COMMON.DBASE'
1673 include 'COMMON.THREAD'
1674 include 'COMMON.TIME1'
1675 dimension itype_pdb(maxres)
1676 common /pizda/ itype_pdb
1678 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1679 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1680 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1682 do i=nstart_sup,nstart_sup+nsup-1
1683 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1684 cd & ' seq_pdb', restyp(itype_pdb(i))
1685 do j=i+2,nstart_sup+nsup-1
1687 ihpb(nhpb)=i+nstart_seq-nstart_sup
1688 jhpb(nhpb)=j+nstart_seq-nstart_sup
1690 dhpb(nhpb)=dist(i,j)
1693 cd write (iout,'(a)') 'Distance constraints:'
1698 cd if (ii.gt.nres) then
1703 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1704 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1705 cd & dhpb(i),forcon(i)
1709 c----------------------------------------------------------------------------
1711 implicit real*8 (a-h,o-z)
1712 include 'DIMENSIONS'
1713 include 'COMMON.MAP'
1714 include 'COMMON.IOUNITS'
1715 character*3 angid(4) /'THE','PHI','ALP','OME'/
1716 character*80 mapcard,ucase
1718 read (inp,'(a)') mapcard
1719 mapcard=ucase(mapcard)
1720 if (index(mapcard,'PHI').gt.0) then
1722 else if (index(mapcard,'THE').gt.0) then
1724 else if (index(mapcard,'ALP').gt.0) then
1726 else if (index(mapcard,'OME').gt.0) then
1729 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1730 stop 'Error - illegal variable spec in MAP card.'
1732 call readi (mapcard,'RES1',res1(imap),0)
1733 call readi (mapcard,'RES2',res2(imap),0)
1734 if (res1(imap).eq.0) then
1735 res1(imap)=res2(imap)
1736 else if (res2(imap).eq.0) then
1737 res2(imap)=res1(imap)
1739 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1741 & 'Error - illegal definition of variable group in MAP.'
1742 stop 'Error - illegal definition of variable group in MAP.'
1744 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1745 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1746 call readi(mapcard,'NSTEP',nstep(imap),0)
1747 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1749 & 'Illegal boundary and/or step size specification in MAP.'
1750 stop 'Illegal boundary and/or step size specification in MAP.'
1755 c----------------------------------------------------------------------------
1756 csa subroutine csaread
1757 csa implicit real*8 (a-h,o-z)
1758 csa include 'DIMENSIONS'
1759 csa include 'COMMON.IOUNITS'
1760 csa include 'COMMON.GEO'
1761 csa include 'COMMON.CSA'
1762 csa include 'COMMON.BANK'
1763 csa include 'COMMON.CONTROL'
1764 csa character*80 ucase
1765 csa character*620 mcmcard
1766 csa call card_concat(mcmcard)
1768 csa call readi(mcmcard,'NCONF',nconf,50)
1769 csa call readi(mcmcard,'NADD',nadd,0)
1770 csa call readi(mcmcard,'JSTART',jstart,1)
1771 csa call readi(mcmcard,'JEND',jend,1)
1772 csa call readi(mcmcard,'NSTMAX',nstmax,500000)
1773 csa call readi(mcmcard,'N0',n0,1)
1774 csa call readi(mcmcard,'N1',n1,6)
1775 csa call readi(mcmcard,'N2',n2,4)
1776 csa call readi(mcmcard,'N3',n3,0)
1777 csa call readi(mcmcard,'N4',n4,0)
1778 csa call readi(mcmcard,'N5',n5,0)
1779 csa call readi(mcmcard,'N6',n6,10)
1780 csa call readi(mcmcard,'N7',n7,0)
1781 csa call readi(mcmcard,'N8',n8,0)
1782 csa call readi(mcmcard,'N9',n9,0)
1783 csa call readi(mcmcard,'N14',n14,0)
1784 csa call readi(mcmcard,'N15',n15,0)
1785 csa call readi(mcmcard,'N16',n16,0)
1786 csa call readi(mcmcard,'N17',n17,0)
1787 csa call readi(mcmcard,'N18',n18,0)
1789 csa vdisulf=(index(mcmcard,'DYNSS').gt.0)
1791 csa call readi(mcmcard,'NDIFF',ndiff,2)
1792 csa call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1793 csa call readi(mcmcard,'IS1',is1,1)
1794 csa call readi(mcmcard,'IS2',is2,8)
1795 csa call readi(mcmcard,'NRAN0',nran0,4)
1796 csa call readi(mcmcard,'NRAN1',nran1,2)
1797 csa call readi(mcmcard,'IRR',irr,1)
1798 csa call readi(mcmcard,'NSEED',nseed,20)
1799 csa call readi(mcmcard,'NTOTAL',ntotal,10000)
1800 csa call reada(mcmcard,'CUT1',cut1,2.0d0)
1801 csa call reada(mcmcard,'CUT2',cut2,5.0d0)
1802 csa call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1803 csa call readi(mcmcard,'ICMAX',icmax,3)
1804 csa call readi(mcmcard,'IRESTART',irestart,0)
1805 csac!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1808 csa call reada(mcmcard,'DELE',dele,20.0d0)
1809 csa call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1810 csa call readi(mcmcard,'IREF',iref,0)
1811 csa call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1812 csa call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1813 csa call readi(mcmcard,'NCONF_IN',nconf_in,0)
1814 csa call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1815 csa write (iout,*) "NCONF_IN",nconf_in
1818 c----------------------------------------------------------------------------
1819 cfmc subroutine mcmfread
1820 cfmc implicit real*8 (a-h,o-z)
1821 cfmc include 'DIMENSIONS'
1822 cfmc include 'COMMON.MCMF'
1823 cfmc include 'COMMON.IOUNITS'
1824 cfmc include 'COMMON.GEO'
1825 cfmc character*80 ucase
1826 cfmc character*620 mcmcard
1827 cfmc call card_concat(mcmcard)
1829 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1830 cfmc write(iout,*)'MAXRANT=',maxrant
1831 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1832 cfmc write(iout,*)'MAXFAM=',maxfam
1833 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1834 cfmc write(iout,*)'NNET1=',nnet1
1835 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1836 cfmc write(iout,*)'NNET2=',nnet2
1837 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1838 cfmc write(iout,*)'NNET3=',nnet3
1839 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1840 cfmc write(iout,*)'ILASTT=',ilastt
1841 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1842 cfmc write(iout,*)'MAXSTR=',maxstr
1843 cfmc maxstr_f=maxstr/maxfam
1844 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1845 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1846 cfmc write(iout,*)'NMCMF=',nmcmf
1847 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1848 cfmc write(iout,*)'IFOCUS=',ifocus
1849 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1850 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1851 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1852 cfmc write(iout,*)'INTPRT=',intprt
1853 cfmc call readi(mcmcard,'IPRT',iprt,100)
1854 cfmc write(iout,*)'IPRT=',iprt
1855 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1856 cfmc write(iout,*)'IMAXTR=',imaxtr
1857 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1858 cfmc write(iout,*)'MAXEVEN=',maxeven
1859 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1860 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1861 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1862 cfmc write(iout,*)'INIMIN=',inimin
1863 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1864 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1865 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1866 cfmc write(iout,*)'NTHREAD=',nthread
1867 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1868 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1869 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1870 cfmc write(iout,*)'MAXPERT=',maxpert
1871 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1872 cfmc write(iout,*)'IRMSD=',irmsd
1873 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1874 cfmc write(iout,*)'DENEMIN=',denemin
1875 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1876 cfmc write(iout,*)'RCUT1S=',rcut1s
1877 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1878 cfmc write(iout,*)'RCUT1E=',rcut1e
1879 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1880 cfmc write(iout,*)'RCUT2S=',rcut2s
1881 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1882 cfmc write(iout,*)'RCUT2E=',rcut2e
1883 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1884 cfmc write(iout,*)'DPERT1=',d_pert1
1885 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1886 cfmc write(iout,*)'DPERT1A=',d_pert1a
1887 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1888 cfmc write(iout,*)'DPERT2=',d_pert2
1889 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1890 cfmc write(iout,*)'DPERT2A=',d_pert2a
1891 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1892 cfmc write(iout,*)'DPERT2B=',d_pert2b
1893 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1894 cfmc write(iout,*)'DPERT2C=',d_pert2c
1895 cfmc d_pert1=deg2rad*d_pert1
1896 cfmc d_pert1a=deg2rad*d_pert1a
1897 cfmc d_pert2=deg2rad*d_pert2
1898 cfmc d_pert2a=deg2rad*d_pert2a
1899 cfmc d_pert2b=deg2rad*d_pert2b
1900 cfmc d_pert2c=deg2rad*d_pert2c
1901 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1902 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1903 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1904 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1905 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1906 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1907 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1908 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1909 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1910 cfmc write(iout,*)'RCUTINI=',rcutini
1911 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1912 cfmc write(iout,*)'GRAT=',grat
1913 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1914 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1918 c----------------------------------------------------------------------------
1920 implicit real*8 (a-h,o-z)
1921 include 'DIMENSIONS'
1922 include 'COMMON.MCM'
1923 include 'COMMON.MCE'
1924 include 'COMMON.IOUNITS'
1926 character*320 mcmcard
1927 call card_concat(mcmcard)
1928 call readi(mcmcard,'MAXACC',maxacc,100)
1929 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1930 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1931 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1932 call readi(mcmcard,'MAXREPM',maxrepm,200)
1933 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1934 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1935 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1936 call reada(mcmcard,'E_UP',e_up,5.0D0)
1937 call reada(mcmcard,'DELTE',delte,0.1D0)
1938 call readi(mcmcard,'NSWEEP',nsweep,5)
1939 call readi(mcmcard,'NSTEPH',nsteph,0)
1940 call readi(mcmcard,'NSTEPC',nstepc,0)
1941 call reada(mcmcard,'TMIN',tmin,298.0D0)
1942 call reada(mcmcard,'TMAX',tmax,298.0D0)
1943 call readi(mcmcard,'NWINDOW',nwindow,0)
1944 call readi(mcmcard,'PRINT_MC',print_mc,0)
1945 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1946 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1947 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1948 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1949 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1950 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1951 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1952 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1953 if (nwindow.gt.0) then
1954 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1956 winlen(i)=winend(i)-winstart(i)+1
1959 if (tmax.lt.tmin) tmax=tmin
1960 if (tmax.eq.tmin) then
1964 if (nstepc.gt.0 .and. nsteph.gt.0) then
1965 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1966 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1968 C Probabilities of different move types
1969 sumpro_type(0)=0.0D0
1970 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1971 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1972 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1973 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1974 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1975 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1976 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1978 print *,'i',i,' sumprotype',sumpro_type(i)
1979 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1980 print *,'i',i,' sumprotype',sumpro_type(i)
1984 c----------------------------------------------------------------------------
1985 subroutine read_minim
1986 implicit real*8 (a-h,o-z)
1987 include 'DIMENSIONS'
1988 include 'COMMON.MINIM'
1989 include 'COMMON.IOUNITS'
1991 character*320 minimcard
1992 call card_concat(minimcard)
1993 call readi(minimcard,'MAXMIN',maxmin,2000)
1994 call readi(minimcard,'MAXFUN',maxfun,5000)
1995 call readi(minimcard,'MINMIN',minmin,maxmin)
1996 call readi(minimcard,'MINFUN',minfun,maxmin)
1997 call reada(minimcard,'TOLF',tolf,1.0D-2)
1998 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1999 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
2000 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
2001 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
2002 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2003 & 'Options in energy minimization:'
2004 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
2005 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
2006 & 'MinMin:',MinMin,' MinFun:',MinFun,
2007 & ' TolF:',TolF,' RTolF:',RTolF
2010 c----------------------------------------------------------------------------
2011 subroutine read_angles(kanal,*)
2012 implicit real*8 (a-h,o-z)
2013 include 'DIMENSIONS'
2014 include 'COMMON.GEO'
2015 include 'COMMON.VAR'
2016 include 'COMMON.CHAIN'
2017 include 'COMMON.IOUNITS'
2018 include 'COMMON.CONTROL'
2019 c Read angles from input
2021 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
2022 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
2023 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
2024 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
2027 c 9/7/01 avoid 180 deg valence angle
2028 if (theta(i).gt.179.99d0) theta(i)=179.99d0
2030 theta(i)=deg2rad*theta(i)
2031 phi(i)=deg2rad*phi(i)
2032 alph(i)=deg2rad*alph(i)
2033 omeg(i)=deg2rad*omeg(i)
2038 c----------------------------------------------------------------------------
2039 subroutine reada(rekord,lancuch,wartosc,default)
2041 character*(*) rekord,lancuch
2042 double precision wartosc,default
2045 iread=index(rekord,lancuch)
2046 if (iread.eq.0) then
2050 iread=iread+ilen(lancuch)+1
2051 read (rekord(iread:),*,err=10,end=10) wartosc
2056 c----------------------------------------------------------------------------
2057 subroutine readi(rekord,lancuch,wartosc,default)
2059 character*(*) rekord,lancuch
2060 integer wartosc,default
2063 iread=index(rekord,lancuch)
2064 if (iread.eq.0) then
2068 iread=iread+ilen(lancuch)+1
2069 read (rekord(iread:),*,err=10,end=10) wartosc
2074 c----------------------------------------------------------------------------
2075 subroutine multreadi(rekord,lancuch,tablica,dim,default)
2078 integer tablica(dim),default
2079 character*(*) rekord,lancuch
2086 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2087 if (iread.eq.0) return
2088 iread=iread+ilen(lancuch)+1
2089 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2092 c----------------------------------------------------------------------------
2093 subroutine multreada(rekord,lancuch,tablica,dim,default)
2096 double precision tablica(dim),default
2097 character*(*) rekord,lancuch
2104 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2105 if (iread.eq.0) return
2106 iread=iread+ilen(lancuch)+1
2107 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2110 c----------------------------------------------------------------------------
2111 subroutine openunits
2112 implicit real*8 (a-h,o-z)
2113 include 'DIMENSIONS'
2116 character*16 form,nodename
2119 include 'COMMON.SETUP'
2120 include 'COMMON.IOUNITS'
2122 include 'COMMON.CONTROL'
2123 integer lenpre,lenpot,ilen,lentmp
2125 character*3 out1file_text,ucase
2128 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2129 call getenv_loc("PREFIX",prefix)
2131 call getenv_loc("POT",pot)
2132 call getenv_loc("DIRTMP",tmpdir)
2133 call getenv_loc("CURDIR",curdir)
2134 call getenv_loc("OUT1FILE",out1file_text)
2135 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2136 out1file_text=ucase(out1file_text)
2137 if (out1file_text(1:1).eq."Y") then
2140 out1file=fg_rank.gt.0
2145 if (lentmp.gt.0) then
2146 write (*,'(80(1h!))')
2147 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2148 write (*,'(80(1h!))')
2149 write (*,*)"All output files will be on node /tmp directory."
2151 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2152 if (me.eq.king) then
2153 write (*,*) "The master node is ",nodename
2154 else if (fg_rank.eq.0) then
2155 write (*,*) "I am the CG slave node ",nodename
2157 write (*,*) "I am the FG slave node ",nodename
2160 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2161 lenpre = lentmp+lenpre+1
2163 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2164 C Get the names and open the input files
2165 #if defined(WINIFL) || defined(WINPGI)
2166 open(1,file=pref_orig(:ilen(pref_orig))//
2167 & '.inp',status='old',readonly,shared)
2168 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2169 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2170 C Get parameter filenames and open the parameter files.
2171 call getenv_loc('BONDPAR',bondname)
2172 open (ibond,file=bondname,status='old',readonly,shared)
2173 call getenv_loc('THETPAR',thetname)
2174 open (ithep,file=thetname,status='old',readonly,shared)
2176 call getenv_loc('THETPARPDB',thetname_pdb)
2177 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2179 call getenv_loc('ROTPAR',rotname)
2180 open (irotam,file=rotname,status='old',readonly,shared)
2182 call getenv_loc('ROTPARPDB',rotname_pdb)
2183 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2185 call getenv_loc('TORPAR',torname)
2186 open (itorp,file=torname,status='old',readonly,shared)
2187 call getenv_loc('TORDPAR',tordname)
2188 open (itordp,file=tordname,status='old',readonly,shared)
2189 call getenv_loc('FOURIER',fouriername)
2190 open (ifourier,file=fouriername,status='old',readonly,shared)
2191 call getenv_loc('ELEPAR',elename)
2192 open (ielep,file=elename,status='old',readonly,shared)
2193 call getenv_loc('SIDEPAR',sidename)
2194 open (isidep,file=sidename,status='old',readonly,shared)
2195 #elif (defined CRAY) || (defined AIX)
2196 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2198 c print *,"Processor",myrank," opened file 1"
2199 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2200 c print *,"Processor",myrank," opened file 9"
2201 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2202 C Get parameter filenames and open the parameter files.
2203 call getenv_loc('BONDPAR',bondname)
2204 open (ibond,file=bondname,status='old',action='read')
2205 c print *,"Processor",myrank," opened file IBOND"
2206 call getenv_loc('THETPAR',thetname)
2207 open (ithep,file=thetname,status='old',action='read')
2208 c print *,"Processor",myrank," opened file ITHEP"
2210 call getenv_loc('THETPARPDB',thetname_pdb)
2211 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2213 call getenv_loc('ROTPAR',rotname)
2214 open (irotam,file=rotname,status='old',action='read')
2215 c print *,"Processor",myrank," opened file IROTAM"
2217 call getenv_loc('ROTPARPDB',rotname_pdb)
2218 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2220 call getenv_loc('TORPAR',torname)
2221 open (itorp,file=torname,status='old',action='read')
2222 c print *,"Processor",myrank," opened file ITORP"
2223 call getenv_loc('TORDPAR',tordname)
2224 open (itordp,file=tordname,status='old',action='read')
2225 c print *,"Processor",myrank," opened file ITORDP"
2226 call getenv_loc('SCCORPAR',sccorname)
2227 open (isccor,file=sccorname,status='old',action='read')
2228 c print *,"Processor",myrank," opened file ISCCOR"
2229 call getenv_loc('FOURIER',fouriername)
2230 open (ifourier,file=fouriername,status='old',action='read')
2231 c print *,"Processor",myrank," opened file IFOURIER"
2232 call getenv_loc('ELEPAR',elename)
2233 open (ielep,file=elename,status='old',action='read')
2234 c print *,"Processor",myrank," opened file IELEP"
2235 call getenv_loc('SIDEPAR',sidename)
2236 open (isidep,file=sidename,status='old',action='read')
2237 c print *,"Processor",myrank," opened file ISIDEP"
2238 c print *,"Processor",myrank," opened parameter files"
2240 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2241 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2242 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2243 C Get parameter filenames and open the parameter files.
2244 call getenv_loc('BONDPAR',bondname)
2245 open (ibond,file=bondname,status='old')
2246 call getenv_loc('THETPAR',thetname)
2247 open (ithep,file=thetname,status='old')
2249 call getenv_loc('THETPARPDB',thetname_pdb)
2250 open (ithep_pdb,file=thetname_pdb,status='old')
2252 call getenv_loc('ROTPAR',rotname)
2253 open (irotam,file=rotname,status='old')
2255 call getenv_loc('ROTPARPDB',rotname_pdb)
2256 open (irotam_pdb,file=rotname_pdb,status='old')
2258 call getenv_loc('TORPAR',torname)
2259 open (itorp,file=torname,status='old')
2260 call getenv_loc('TORDPAR',tordname)
2261 open (itordp,file=tordname,status='old')
2262 call getenv_loc('SCCORPAR',sccorname)
2263 open (isccor,file=sccorname,status='old')
2264 call getenv_loc('FOURIER',fouriername)
2265 open (ifourier,file=fouriername,status='old')
2266 call getenv_loc('ELEPAR',elename)
2267 open (ielep,file=elename,status='old')
2268 call getenv_loc('SIDEPAR',sidename)
2269 open (isidep,file=sidename,status='old')
2271 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2273 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2274 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2275 C Get parameter filenames and open the parameter files.
2276 call getenv_loc('BONDPAR',bondname)
2277 open (ibond,file=bondname,status='old',action='read')
2278 call getenv_loc('THETPAR',thetname)
2279 open (ithep,file=thetname,status='old',action='read')
2281 call getenv_loc('THETPARPDB',thetname_pdb)
2282 print *,"thetname_pdb ",thetname_pdb
2283 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2284 print *,ithep_pdb," opened"
2286 call getenv_loc('ROTPAR',rotname)
2287 open (irotam,file=rotname,status='old',action='read')
2289 call getenv_loc('ROTPARPDB',rotname_pdb)
2290 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2292 call getenv_loc('TORPAR',torname)
2293 open (itorp,file=torname,status='old',action='read')
2294 call getenv_loc('TORDPAR',tordname)
2295 open (itordp,file=tordname,status='old',action='read')
2296 call getenv_loc('SCCORPAR',sccorname)
2297 open (isccor,file=sccorname,status='old',action='read')
2298 call getenv_loc('FOURIER',fouriername)
2299 open (ifourier,file=fouriername,status='old',action='read')
2300 call getenv_loc('ELEPAR',elename)
2301 open (ielep,file=elename,status='old',action='read')
2302 call getenv_loc('SIDEPAR',sidename)
2303 open (isidep,file=sidename,status='old',action='read')
2307 C 8/9/01 In the newest version SCp interaction constants are read from a file
2308 C Use -DOLDSCP to use hard-coded constants instead.
2310 call getenv_loc('SCPPAR',scpname)
2311 #if defined(WINIFL) || defined(WINPGI)
2312 open (iscpp,file=scpname,status='old',readonly,shared)
2313 #elif (defined CRAY) || (defined AIX)
2314 open (iscpp,file=scpname,status='old',action='read')
2316 open (iscpp,file=scpname,status='old')
2318 open (iscpp,file=scpname,status='old',action='read')
2321 call getenv_loc('PATTERN',patname)
2322 #if defined(WINIFL) || defined(WINPGI)
2323 open (icbase,file=patname,status='old',readonly,shared)
2324 #elif (defined CRAY) || (defined AIX)
2325 open (icbase,file=patname,status='old',action='read')
2327 open (icbase,file=patname,status='old')
2329 open (icbase,file=patname,status='old',action='read')
2332 C Open output file only for CG processes
2333 c print *,"Processor",myrank," fg_rank",fg_rank
2334 if (fg_rank.eq.0) then
2336 if (nodes.eq.1) then
2339 npos = dlog10(dfloat(nodes-1))+1
2341 if (npos.lt.3) npos=3
2342 write (liczba,'(i1)') npos
2343 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2345 write (liczba,form) me
2346 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2347 & liczba(:ilen(liczba))
2348 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2350 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2352 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2353 & liczba(:ilen(liczba))//'.mol2'
2354 cartname=prefix(:lenpre)//'_'//pot(:lenpot)//
2355 & liczba(:ilen(liczba))//'.x'
2356 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2357 & liczba(:ilen(liczba))//'.stat'
2359 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2360 & //liczba(:ilen(liczba))//'.stat')
2361 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2364 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2365 & liczba(:ilen(liczba))//'.const'
2370 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2371 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2372 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2373 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2374 cartname=prefix(:lenpre)//'_'//pot(:lenpot)//'.x'
2375 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2377 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2379 rest2name=prefix(:ilen(prefix))//'.rst'
2381 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2384 #if defined(AIX) || defined(PGI)
2385 if (me.eq.king .or. .not. out1file)
2386 & open(iout,file=outname,status='unknown')
2389 if (fg_rank.gt.0) then
2390 write (liczba,'(i3.3)') myrank/nfgtasks
2391 write (ll,'(bz,i3.3)') fg_rank
2392 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2398 open(igeom,file=intname,status='unknown',position='append')
2399 open(ipdb,file=pdbname,status='unknown')
2400 open(imol2,file=mol2name,status='unknown')
2401 open(istat,file=statname,status='unknown',position='append')
2403 c1out open(iout,file=outname,status='unknown')
2406 if (me.eq.king .or. .not.out1file)
2407 & open(iout,file=outname,status='unknown')
2410 if (fg_rank.gt.0) then
2411 print "Processor",fg_rank," opening output file"
2412 write (liczba,'(i3.3)') myrank/nfgtasks
2413 write (ll,'(bz,i3.3)') fg_rank
2414 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2420 open(igeom,file=intname,status='unknown',access='append')
2421 open(ipdb,file=pdbname,status='unknown')
2422 open(imol2,file=mol2name,status='unknown')
2423 open(istat,file=statname,status='unknown',access='append')
2425 c1out open(iout,file=outname,status='unknown')
2428 csa csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2429 csa csa_seed=prefix(:lenpre)//'.CSA.seed'
2430 csa csa_history=prefix(:lenpre)//'.CSA.history'
2431 csa csa_bank=prefix(:lenpre)//'.CSA.bank'
2432 csa csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2433 csa csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2434 csa csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2435 csac!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2436 csa csa_int=prefix(:lenpre)//'.int'
2437 csa csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2438 csa csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2439 csa csa_in=prefix(:lenpre)//'.CSA.in'
2440 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2443 write (iout,'(80(1h-))')
2444 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2445 write (iout,'(80(1h-))')
2446 write (iout,*) "Input file : ",
2447 & pref_orig(:ilen(pref_orig))//'.inp'
2448 write (iout,*) "Output file : ",
2449 & outname(:ilen(outname))
2451 write (iout,*) "Sidechain potential file : ",
2452 & sidename(:ilen(sidename))
2454 write (iout,*) "SCp potential file : ",
2455 & scpname(:ilen(scpname))
2457 write (iout,*) "Electrostatic potential file : ",
2458 & elename(:ilen(elename))
2459 write (iout,*) "Cumulant coefficient file : ",
2460 & fouriername(:ilen(fouriername))
2461 write (iout,*) "Torsional parameter file : ",
2462 & torname(:ilen(torname))
2463 write (iout,*) "Double torsional parameter file : ",
2464 & tordname(:ilen(tordname))
2465 write (iout,*) "SCCOR parameter file : ",
2466 & sccorname(:ilen(sccorname))
2467 write (iout,*) "Bond & inertia constant file : ",
2468 & bondname(:ilen(bondname))
2469 write (iout,*) "Bending parameter file : ",
2470 & thetname(:ilen(thetname))
2471 write (iout,*) "Rotamer parameter file : ",
2472 & rotname(:ilen(rotname))
2473 write (iout,*) "Threading database : ",
2474 & patname(:ilen(patname))
2476 &write (iout,*)" DIRTMP : ",
2478 write (iout,'(80(1h-))')
2482 c----------------------------------------------------------------------------
2483 subroutine card_concat(card)
2484 implicit real*8 (a-h,o-z)
2485 include 'DIMENSIONS'
2486 include 'COMMON.IOUNITS'
2488 character*80 karta,ucase
2490 read (inp,'(a)') karta
2493 do while (karta(80:80).eq.'&')
2494 card=card(:ilen(card)+1)//karta(:79)
2495 read (inp,'(a)') karta
2498 card=card(:ilen(card)+1)//karta
2501 c----------------------------------------------------------------------------------
2503 implicit real*8 (a-h,o-z)
2504 include 'DIMENSIONS'
2505 include 'COMMON.CHAIN'
2506 include 'COMMON.IOUNITS'
2508 include 'COMMON.CONTROL'
2509 open(irest2,file=rest2name,status='unknown')
2510 read(irest2,*) totT,EK,potE,totE,t_bath
2512 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2515 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2517 if(usampl.or.homol_nset.gt.1) then
2518 read (irest2,*) iset
2523 c---------------------------------------------------------------------------------
2524 subroutine read_fragments
2525 implicit real*8 (a-h,o-z)
2526 include 'DIMENSIONS'
2530 include 'COMMON.SETUP'
2531 include 'COMMON.CHAIN'
2532 include 'COMMON.IOUNITS'
2534 include 'COMMON.CONTROL'
2536 read(inp,*) nset,nfrag,npair,nfrag_back
2537 if(me.eq.king.or..not.out1file)
2538 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2539 & " nfrag_back",nfrag_back
2541 read(inp,*) mset(iset1)
2543 read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
2545 if(me.eq.king.or..not.out1file)
2546 & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
2547 & ifrag(2,i,iset1), qinfrag(i,iset1)
2550 read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
2552 if(me.eq.king.or..not.out1file)
2553 & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
2554 & ipair(2,i,iset1), qinpair(i,iset1)
2557 read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
2558 & wfrag_back(3,i,iset1),
2559 & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
2560 if(me.eq.king.or..not.out1file)
2561 & write(iout,*) "A",i,wfrag_back(1,i,iset1),
2562 & wfrag_back(2,i,iset1),
2563 & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
2564 & ifrag_back(2,i,iset1)
2569 c-------------------------------------------------------------------------------
2570 subroutine read_dist_constr
2571 implicit real*8 (a-h,o-z)
2572 include 'DIMENSIONS'
2576 include 'COMMON.SETUP'
2577 include 'COMMON.CONTROL'
2578 include 'COMMON.CHAIN'
2579 include 'COMMON.IOUNITS'
2580 include 'COMMON.SBRIDGE'
2581 integer ifrag_(2,100),ipair_(2,100)
2582 double precision wfrag_(100),wpair_(100)
2583 character*500 controlcard
2584 logical normalize,next
2586 double precision xlink(4,0:4) /
2588 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
2589 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
2590 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
2591 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
2592 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
2593 c print *, "WCHODZE"
2594 write (iout,*) "Calling read_dist_constr"
2595 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2601 call card_concat(controlcard)
2602 next = index(controlcard,"NEXT").gt.0
2603 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2604 write (iout,*) "restr_type",restr_type
2605 call readi(controlcard,"NFRAG",nfrag_,0)
2606 call readi(controlcard,"NFRAG",nfrag_,0)
2607 call readi(controlcard,"NPAIR",npair_,0)
2608 call readi(controlcard,"NDIST",ndist_,0)
2609 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2610 if (restr_type.eq.10)
2611 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2612 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2613 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2614 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2615 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2616 normalize = index(controlcard,"NORMALIZE").gt.0
2617 write (iout,*) "WBOLTZD",wboltzd
2618 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2619 c write (iout,*) "IFRAG"
2621 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2623 c write (iout,*) "IPAIR"
2625 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2627 if (nfrag_.gt.0) write (iout,*)
2628 & "Distance restraints as generated from reference structure"
2630 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2631 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2632 & ifrag_(2,i)=nstart_sup+nsup-1
2633 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2635 if (wfrag_(i).eq.0.0d0) cycle
2636 do j=ifrag_(1,i),ifrag_(2,i)-1
2637 do k=j+1,ifrag_(2,i)
2638 c write (iout,*) "j",j," k",k
2640 if (restr_type.eq.1) then
2646 forcon(nhpb)=wfrag_(i)
2647 else if (constr_dist.eq.2) then
2648 if (ddjk.le.dist_cut) then
2654 forcon(nhpb)=wfrag_(i)
2656 else if (restr_type.eq.3) then
2662 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2665 if (.not.out1file .or. me.eq.king)
2666 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2667 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2669 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2670 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2676 if (wpair_(i).eq.0.0d0) cycle
2684 do j=ifrag_(1,ii),ifrag_(2,ii)
2685 do k=ifrag_(1,jj),ifrag_(2,jj)
2686 if (restr_type.eq.1) then
2692 forcon(nhpb)=wfrag_(i)
2693 else if (constr_dist.eq.2) then
2694 if (ddjk.le.dist_cut) then
2700 forcon(nhpb)=wfrag_(i)
2702 else if (restr_type.eq.3) then
2708 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2711 if (.not.out1file .or. me.eq.king)
2712 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
2713 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2715 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.restr ",
2716 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2723 write (iout,*) "Distance restraints as read from input"
2725 if (restr_type.eq.11) then
2726 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2727 & dhpb1(nhpb+1),ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2728 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2729 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2731 irestr_type(nhpb)=11
2733 if (.not.out1file .or. me.eq.king)
2734 & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2735 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2736 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2738 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2739 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2740 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2742 if (ibecarb(nhpb).gt.0) then
2743 ihpb(nhpb)=ihpb(nhpb)+nres
2744 jhpb(nhpb)=jhpb(nhpb)+nres
2746 else if (constr_dist.eq.10) then
2747 c Cross-lonk Markov-like potential
2748 call card_concat(controlcard)
2749 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2750 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2752 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2753 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2754 if (index(controlcard,"ZL").gt.0) then
2756 else if (index(controlcard,"ADH").gt.0) then
2758 else if (index(controlcard,"PDH").gt.0) then
2760 else if (index(controlcard,"DSS").gt.0) then
2765 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2766 & xlink(1,link_type))
2767 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2768 & xlink(2,link_type))
2769 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2770 & xlink(3,link_type))
2771 call reada(controlcard,"SIGMA",forcon(nhpb+1),
2772 & xlink(4,link_type))
2773 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2774 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2775 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2776 if (forcon(nhpb+1).le.0.0d0 .or.
2777 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2779 irestr_type(nhpb)=10
2780 if (ibecarb(nhpb).gt.0) then
2781 ihpb(nhpb)=ihpb(nhpb)+nres
2782 jhpb(nhpb)=jhpb(nhpb)+nres
2785 if (.not.out1file .or. me.eq.king)
2786 & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2787 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2788 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2791 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2792 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2793 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2798 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2799 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2800 if (forcon(nhpb+1).gt.0.0d0) then
2802 if (dhpb1(nhpb).eq.0.0d0) then
2807 if (ibecarb(nhpb).gt.0) then
2808 ihpb(nhpb)=ihpb(nhpb)+nres
2809 jhpb(nhpb)=jhpb(nhpb)+nres
2811 if (dhpb(nhpb).eq.0.0d0)
2812 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2815 if (.not.out1file .or. me.eq.king)
2816 & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2817 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2819 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2820 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2823 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2824 C if (forcon(nhpb+1).gt.0.0d0) then
2826 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2834 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
2835 & fordepthmax=fordepth(i)
2838 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
2843 c-------------------------------------------------------------------------------
2845 subroutine read_constr_homology
2847 include 'DIMENSIONS'
2851 include 'COMMON.SETUP'
2852 include 'COMMON.CONTROL'
2853 include 'COMMON.CHAIN'
2854 include 'COMMON.IOUNITS'
2856 include 'COMMON.GEO'
2857 include 'COMMON.INTERACT'
2858 include 'COMMON.NAMES'
2860 c For new homol impl
2862 include 'COMMON.VAR'
2865 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2867 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2868 c & sigma_odl_temp(maxres,maxres,max_template)
2870 character*24 model_ki_dist, model_ki_angle
2871 character*500 controlcard
2872 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
2873 logical lprn /.true./
2878 c FP - Nov. 2014 Temporary specifications for new vars
2880 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
2881 double precision, dimension (max_template,maxres) :: rescore
2882 double precision, dimension (max_template,maxres) :: rescore2
2883 character*24 pdbfile,tpl_k_rescore
2884 c -----------------------------------------------------------------
2885 c Reading multiple PDB ref structures and calculation of retraints
2886 c not using pre-computed ones stored in files model_ki_{dist,angle}
2888 c -----------------------------------------------------------------
2891 c Alternative: reading from input
2892 call card_concat(controlcard)
2893 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2894 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2895 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2896 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2897 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2898 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
2899 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2900 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
2901 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
2902 if(.not.read2sigma.and.start_from_model) then
2903 write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
2904 start_from_model=.false.
2906 if(start_from_model) write(iout,*) 'START_FROM_MODELS is ON'
2907 if(start_from_model .and. rest) then
2908 write(iout,*) 'START_FROM_MODELS is OFF'
2909 write(iout,*) 'remove restart keyword from input'
2911 if (homol_nset.gt.1)then
2912 call card_concat(controlcard)
2913 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2914 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2915 write(iout,*) "iset homology_weight "
2917 write(iout,*) i,waga_homology(i)
2920 iset=mod(kolor,homol_nset)+1
2923 waga_homology(1)=1.0
2926 cd write (iout,*) "nnt",nnt," nct",nct
2933 write(iout,*) 'nnt=',nnt,'nct=',nct
2936 do k=1,constr_homology
2949 if (read_homol_frag) then
2950 call read_klapaucjusz
2953 do k=1,constr_homology
2955 read(inp,'(a)') pdbfile
2956 c Next stament causes error upon compilation (?)
2957 c if(me.eq.king.or. .not. out1file)
2958 c write (iout,'(2a)') 'PDB data will be read from file ',
2959 c & pdbfile(:ilen(pdbfile))
2960 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
2961 & pdbfile(:ilen(pdbfile))
2962 open(ipdbin,file=pdbfile,status='old',err=33)
2964 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
2965 & pdbfile(:ilen(pdbfile))
2968 c print *,'Begin reading pdb data'
2970 c Files containing res sim or local scores (former containing sigmas)
2973 write(kic2,'(bz,i2.2)') k
2975 tpl_k_rescore="template"//kic2//".sco"
2978 if (read2sigma) then
2979 call readpdb_template(k)
2984 c Distance restraints
2987 C Copy the coordinates from reference coordinates (?)
2991 c write (iout,*) "c(",j,i,") =",c(j,i)
2995 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
2997 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
2998 open (ientin,file=tpl_k_rescore,status='old')
2999 if (nnt.gt.1) rescore(k,1)=0.0d0
3000 do irec=nnt,nct ! loop for reading res sim
3001 if (read2sigma) then
3002 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3005 idomain(k,i_tmp)=idomain_tmp
3006 rescore(k,i_tmp)=rescore_tmp
3007 rescore2(k,i_tmp)=rescore2_tmp
3008 write(iout,'(a7,i5,2f10.5,i5)') "rescore",
3009 & i_tmp,rescore2_tmp,rescore_tmp,
3013 read (ientin,*,end=1401) rescore_tmp
3015 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3016 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3017 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3022 if (waga_dist.ne.0.0d0) then
3030 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3031 c write (iout,*) k,i,j,distal,dist2_cut
3033 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3034 & .and. distal.le.dist2_cut ) then
3040 c write (iout,*) "k",k
3041 c write (iout,*) "i",i," j",j," constr_homology",
3046 if (read2sigma) then
3049 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3051 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3052 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3053 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3055 if (odl(k,ii).le.dist_cut) then
3056 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3059 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3060 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3062 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3063 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3067 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3070 l_homo(k,ii)=.false.
3077 c Theta, dihedral and SC retraints
3079 if (waga_angle.gt.0.0d0) then
3080 c open (ientin,file=tpl_k_sigma_dih,status='old')
3081 c do irec=1,maxres-3 ! loop for reading sigma_dih
3082 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3083 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3084 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3085 c & sigma_dih(k,i+nnt-1)
3090 if (idomain(k,i).eq.0) then
3094 dih(k,i)=phiref(i) ! right?
3095 c read (ientin,*) sigma_dih(k,i) ! original variant
3096 c write (iout,*) "dih(",k,i,") =",dih(k,i)
3097 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3098 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3099 c & "rescore(",k,i-2,") =",rescore(k,i-2),
3100 c & "rescore(",k,i-3,") =",rescore(k,i-3)
3102 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3103 & rescore(k,i-2)+rescore(k,i-3))/4.0
3104 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3105 c write (iout,*) "Raw sigmas for dihedral angle restraints"
3106 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3107 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3108 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
3109 c Instead of res sim other local measure of b/b str reliability possible
3110 if (sigma_dih(k,i).ne.0)
3111 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3112 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3117 if (waga_theta.gt.0.0d0) then
3118 c open (ientin,file=tpl_k_sigma_theta,status='old')
3119 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3120 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3121 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3122 c & sigma_theta(k,i+nnt-1)
3127 do i = nnt+2,nct ! right? without parallel.
3128 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3129 c do i=ithet_start,ithet_end ! with FG parallel.
3130 if (idomain(k,i).eq.0) then
3131 sigma_theta(k,i)=0.0
3134 thetatpl(k,i)=thetaref(i)
3135 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3136 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3137 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3138 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3139 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3140 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3141 & rescore(k,i-2))/3.0
3142 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3143 if (sigma_theta(k,i).ne.0)
3144 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3146 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3147 c rescore(k,i-2) ! right expression ?
3148 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3152 if (waga_d.gt.0.0d0) then
3153 c open (ientin,file=tpl_k_sigma_d,status='old')
3154 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3155 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3156 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3157 c & sigma_d(k,i+nnt-1)
3161 do i = nnt,nct ! right? without parallel.
3162 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3163 c do i=loc_start,loc_end ! with FG parallel.
3164 if (itype(i).eq.10) cycle
3165 if (idomain(k,i).eq.0 ) then
3172 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3173 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3174 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3175 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3176 sigma_d(k,i)=rescore(k,i) ! right expression ?
3177 if (sigma_d(k,i).ne.0)
3178 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3180 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3181 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3182 c read (ientin,*) sigma_d(k,i) ! 1st variant
3183 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
3188 c remove distance restraints not used in any model from the list
3189 c shift data in all arrays
3191 if (waga_dist.ne.0.0d0) then
3197 if (ii_in_use(ii).eq.0.and.liiflag) then
3201 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3202 & .not.liiflag.and.ii.eq.lim_odl) then
3203 if (ii.eq.lim_odl) then
3204 iishift=ii-iistart+1
3209 do ki=iistart,lim_odl-iishift
3210 ires_homo(ki)=ires_homo(ki+iishift)
3211 jres_homo(ki)=jres_homo(ki+iishift)
3212 ii_in_use(ki)=ii_in_use(ki+iishift)
3213 do k=1,constr_homology
3214 odl(k,ki)=odl(k,ki+iishift)
3215 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3216 l_homo(k,ki)=l_homo(k,ki+iishift)
3220 lim_odl=lim_odl-iishift
3226 endif ! .not. klapaucjusz
3228 if (constr_homology.gt.0) call homology_partition
3229 if (constr_homology.gt.0) call init_int_table
3230 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
3231 cd & "lim_xx=",lim_xx
3232 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3233 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3237 if (.not.lprn) return
3238 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3239 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3240 write (iout,*) "Distance restraints from templates"
3242 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3243 & ii,ires_homo(ii),jres_homo(ii),
3244 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3245 & ki=1,constr_homology)
3247 write (iout,*) "Dihedral angle restraints from templates"
3249 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3250 & (rad2deg*dih(ki,i),
3251 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3253 write (iout,*) "Virtual-bond angle restraints from templates"
3255 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3256 & (rad2deg*thetatpl(ki,i),
3257 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3259 write (iout,*) "SC restraints from templates"
3261 write(iout,'(i5,100(4f8.2,4x))') i,
3262 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3263 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3266 c -----------------------------------------------------------------
3269 c -----------------------------------------------------------------
3270 subroutine read_klapaucjusz
3272 include 'DIMENSIONS'
3276 include 'COMMON.SETUP'
3277 include 'COMMON.CONTROL'
3278 include 'COMMON.CHAIN'
3279 include 'COMMON.IOUNITS'
3281 include 'COMMON.GEO'
3282 include 'COMMON.INTERACT'
3283 include 'COMMON.NAMES'
3284 character*256 fragfile
3285 integer ninclust(maxclust),inclust(max_template,maxclust),
3286 & nresclust(maxclust),iresclust(maxres,maxclust)
3289 character*24 model_ki_dist, model_ki_angle
3290 character*500 controlcard
3291 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
3292 logical lprn /.true./
3298 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3299 double precision, dimension (max_template,maxres) :: rescore
3300 double precision, dimension (max_template,maxres) :: rescore2
3301 character*24 pdbfile,tpl_k_rescore
3304 c For new homol impl
3306 include 'COMMON.VAR'
3308 call getenv("FRAGFILE",fragfile)
3309 open(ientin,file=fragfile,status="old",err=10)
3310 read(ientin,*) constr_homology,nclust
3316 do k=1,constr_homology
3317 read(ientin,'(a)') pdbfile
3318 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3319 & pdbfile(:ilen(pdbfile))
3320 open(ipdbin,file=pdbfile,status='old',err=33)
3322 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3323 & pdbfile(:ilen(pdbfile))
3327 call readpdb_template(k)
3335 read(ientin,*) ninclust(i),nresclust(i)
3336 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3337 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3340 c Loop over clusters
3343 do ll = 1,ninclust(l)
3351 idomain(k,iresclust(i,l)+1) = 1
3353 idomain(k,iresclust(i,l)) = 1
3357 c Distance restraints
3360 C Copy the coordinates from reference coordinates (?)
3364 c write (iout,*) "c(",j,i,") =",c(j,i)
3367 call int_from_cart(.true.,.false.)
3368 call sc_loc_geom(.false.)
3370 thetaref(i)=theta(i)
3373 if (waga_dist.ne.0.0d0) then
3381 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3382 c write (iout,*) k,i,j,distal,dist2_cut
3384 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3385 & .and. distal.le.dist2_cut ) then
3391 c write (iout,*) "k",k
3392 c write (iout,*) "i",i," j",j," constr_homology",
3397 if (read2sigma) then
3400 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3402 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3403 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3404 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3406 if (odl(k,ii).le.dist_cut) then
3407 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3410 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3411 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3413 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3414 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3418 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3421 c l_homo(k,ii)=.false.
3428 c Theta, dihedral and SC retraints
3430 if (waga_angle.gt.0.0d0) then
3432 if (idomain(k,i).eq.0) then
3433 c sigma_dih(k,i)=0.0
3437 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3438 & rescore(k,i-2)+rescore(k,i-3))/4.0
3439 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3440 c & " sigma_dihed",sigma_dih(k,i)
3441 if (sigma_dih(k,i).ne.0)
3442 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3447 if (waga_theta.gt.0.0d0) then
3449 if (idomain(k,i).eq.0) then
3450 c sigma_theta(k,i)=0.0
3453 thetatpl(k,i)=thetaref(i)
3454 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3455 & rescore(k,i-2))/3.0
3456 if (sigma_theta(k,i).ne.0)
3457 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3461 if (waga_d.gt.0.0d0) then
3463 if (itype(i).eq.10) cycle
3464 if (idomain(k,i).eq.0 ) then
3471 sigma_d(k,i)=rescore(k,i)
3472 if (sigma_d(k,i).ne.0)
3473 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3474 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3480 c remove distance restraints not used in any model from the list
3481 c shift data in all arrays
3483 if (waga_dist.ne.0.0d0) then
3489 if (ii_in_use(ii).eq.0.and.liiflag) then
3493 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3494 & .not.liiflag.and.ii.eq.lim_odl) then
3495 if (ii.eq.lim_odl) then
3496 iishift=ii-iistart+1
3501 do ki=iistart,lim_odl-iishift
3502 ires_homo(ki)=ires_homo(ki+iishift)
3503 jres_homo(ki)=jres_homo(ki+iishift)
3504 ii_in_use(ki)=ii_in_use(ki+iishift)
3505 do k=1,constr_homology
3506 odl(k,ki)=odl(k,ki+iishift)
3507 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3508 l_homo(k,ki)=l_homo(k,ki+iishift)
3512 lim_odl=lim_odl-iishift
3519 10 stop "Error infragment file"
3521 c----------------------------------------------------------------------
3524 subroutine flush(iu)
3529 subroutine flush(iu)
3535 c------------------------------------------------------------------------------
3536 subroutine copy_to_tmp(source)
3537 include "DIMENSIONS"
3538 include "COMMON.IOUNITS"
3539 character*(*) source
3540 character* 256 tmpfile
3544 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3545 inquire(file=tmpfile,exist=ex)
3547 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3548 & " to temporary directory..."
3549 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3550 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3554 c------------------------------------------------------------------------------
3555 subroutine move_from_tmp(source)
3556 include "DIMENSIONS"
3557 include "COMMON.IOUNITS"
3558 character*(*) source
3561 write (*,*) "Moving ",source(:ilen(source)),
3562 & " from temporary directory to working directory"
3563 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3564 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3567 c------------------------------------------------------------------------------
3568 subroutine random_init(seed)
3570 C Initialize random number generator
3572 implicit real*8 (a-h,o-z)
3573 include 'DIMENSIONS'
3579 logical OKRandom, prng_restart
3581 integer iseed_array(4)
3583 include 'COMMON.IOUNITS'
3584 include 'COMMON.TIME1'
3585 include 'COMMON.THREAD'
3586 include 'COMMON.SBRIDGE'
3587 include 'COMMON.CONTROL'
3588 include 'COMMON.MCM'
3589 include 'COMMON.MAP'
3590 include 'COMMON.HEADER'
3591 csa include 'COMMON.CSA'
3592 include 'COMMON.CHAIN'
3593 include 'COMMON.MUCA'
3595 include 'COMMON.FFIELD'
3596 include 'COMMON.SETUP'
3597 iseed=-dint(dabs(seed))
3598 if (iseed.eq.0) then
3599 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3600 & 'Random seed undefined. The program will stop.'
3601 write (*,'(/80(1h*)/20x,a/80(1h*))')
3602 & 'Random seed undefined. The program will stop.'
3604 call mpi_finalize(mpi_comm_world,ierr)
3606 stop 'Bad random seed.'
3609 if (fg_rank.eq.0) then
3613 if(me.eq.king .or. .not. out1file)
3614 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3615 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3616 OKRandom = prng_restart(me,iseedi8)
3619 tmp=65536.0d0**(4-i)
3620 iseed_array(i) = dint(seed/tmp)
3621 seed=seed-iseed_array(i)*tmp
3623 if(me.eq.king .or. .not. out1file)
3624 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3625 & (iseed_array(i),i=1,4)
3626 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3627 & (iseed_array(i),i=1,4)
3628 OKRandom = prng_restart(me,iseed_array)
3631 r1=ran_number(0.0D0,1.0D0)
3632 if(me.eq.king .or. .not. out1file)
3633 & write (iout,*) 'ran_num',r1
3634 if (r1.lt.0.0d0) OKRandom=.false.
3636 if (.not.OKRandom) then
3637 write (iout,*) 'PRNG IS NOT WORKING!!!'
3638 print *,'PRNG IS NOT WORKING!!!'
3641 call mpi_abort(mpi_comm_world,error_msg,ierr)
3644 write (iout,*) 'too many processors for parallel prng'
3645 write (*,*) 'too many processors for parallel prng'
3653 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)