2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.SPLITELE'
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 if (modecalc.eq.8) then
33 inquire (file="fort.40",exist=file_exist)
34 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),
57 write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
58 & "The following",npeak,
59 & " NMR peak restraints have been imposed:",
60 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
61 & " score"," type"," ipeak"
63 do j=ipeak(1,i),ipeak(2,i)
64 write (iout,'(5i5,2f8.2,2f10.5,i5)')i,j,ihpb_peak(j),
65 & jhpb_peak(j),ibecarb_peak(j),dhpb_peak(j),dhpb1_peak(j),
66 & forcon_peak(j),fordepth_peak(i),irestr_type_peak(j)
69 write (iout,*) "The ipeak array"
71 write (iout,'(3i5)' ) i,ipeak(1,i),ipeak(2,i)
77 c print *,"Processor",myrank," leaves READRTNS"
80 C-------------------------------------------------------------------------------
81 subroutine read_control
85 implicit real*8 (a-h,o-z)
89 logical OKRandom, prng_restart
92 include 'COMMON.IOUNITS'
93 include 'COMMON.TIME1'
94 include 'COMMON.THREAD'
95 include 'COMMON.SBRIDGE'
96 include 'COMMON.CONTROL'
99 include 'COMMON.HEADER'
101 include 'COMMON.CHAIN'
102 include 'COMMON.MUCA'
104 include 'COMMON.FFIELD'
105 include 'COMMON.INTERACT'
106 include 'COMMON.SETUP'
107 include 'COMMON.SPLITELE'
108 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
109 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
111 character*320 controlcard
116 read (INP,'(a)') titel
117 call card_concat(controlcard)
118 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
119 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
120 call reada(controlcard,'SEED',seed,0.0D0)
121 call random_init(seed)
122 C Set up the time limit (caution! The time must be input in minutes!)
123 read_cart=index(controlcard,'READ_CART').gt.0
124 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
125 write (iout,*) "constr_dist",constr_dist
126 call readi(controlcard,'NSAXS',nsaxs,0)
127 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
128 call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
129 call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
130 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
131 & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
132 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
133 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
134 call readi(controlcard,'SYM',symetr,1)
135 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
136 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
137 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
138 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
139 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
140 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
141 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
142 call reada(controlcard,'DRMS',drms,0.1D0)
143 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
144 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
145 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
146 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
147 write (iout,'(a,f10.1)')'DRMS = ',drms
148 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
149 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
151 call readi(controlcard,'NZ_START',nz_start,0)
152 call readi(controlcard,'NZ_END',nz_end,0)
153 call readi(controlcard,'IZ_SC',iz_sc,0)
155 safety = 60.0d0*safety
158 call reada(controlcard,"T_BATH",t_bath,300.0d0)
159 minim=(index(controlcard,'MINIMIZE').gt.0)
160 dccart=(index(controlcard,'CART').gt.0)
161 overlapsc=(index(controlcard,'OVERLAP').gt.0)
162 overlapsc=.not.overlapsc
163 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
164 searchsc=.not.searchsc
165 sideadd=(index(controlcard,'SIDEADD').gt.0)
166 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
167 outpdb=(index(controlcard,'PDBOUT').gt.0)
168 outmol2=(index(controlcard,'MOL2OUT').gt.0)
169 pdbref=(index(controlcard,'PDBREF').gt.0)
170 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
171 indpdb=index(controlcard,'PDBSTART')
172 extconf=(index(controlcard,'EXTCONF').gt.0)
173 AFMlog=(index(controlcard,'AFM'))
174 selfguide=(index(controlcard,'SELFGUIDE'))
175 c print *,'AFMlog',AFMlog,selfguide,"KUPA"
176 call readi(controlcard,'IPRINT',iprint,0)
177 call readi(controlcard,'MAXGEN',maxgen,10000)
178 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
179 call readi(controlcard,"KDIAG",kdiag,0)
180 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
181 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
182 & write (iout,*) "RESCALE_MODE",rescale_mode
183 split_ene=index(controlcard,'SPLIT_ENE').gt.0
184 if (index(controlcard,'REGULAR').gt.0.0D0) then
185 call reada(controlcard,'WEIDIS',weidis,0.1D0)
189 if (index(controlcard,'CHECKGRAD').gt.0) then
191 if (index(controlcard,'CART').gt.0) then
193 elseif (index(controlcard,'CARINT').gt.0) then
198 elseif (index(controlcard,'THREAD').gt.0) then
200 call readi(controlcard,'THREAD',nthread,0)
201 if (nthread.gt.0) then
202 call reada(controlcard,'WEIDIS',weidis,0.1D0)
205 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
206 stop 'Error termination in Read_Control.'
208 else if (index(controlcard,'MCMA').gt.0) then
210 else if (index(controlcard,'MCEE').gt.0) then
212 else if (index(controlcard,'MULTCONF').gt.0) then
214 else if (index(controlcard,'MAP').gt.0) then
216 call readi(controlcard,'MAP',nmap,0)
217 else if (index(controlcard,'CSA').gt.0) then
219 crc else if (index(controlcard,'ZSCORE').gt.0) then
221 crc ZSCORE is rm from UNRES, modecalc=9 is available
224 cfcm else if (index(controlcard,'MCMF').gt.0) then
226 else if (index(controlcard,'SOFTREG').gt.0) then
228 else if (index(controlcard,'CHECK_BOND').gt.0) then
230 else if (index(controlcard,'TEST').gt.0) then
232 else if (index(controlcard,'MD').gt.0) then
234 else if (index(controlcard,'RE ').gt.0) then
238 lmuca=index(controlcard,'MUCA').gt.0
239 call readi(controlcard,'MUCADYN',mucadyn,0)
240 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
241 if (lmuca .and. (me.eq.king .or. .not.out1file ))
243 write (iout,*) 'MUCADYN=',mucadyn
244 write (iout,*) 'MUCASMOOTH=',muca_smooth
247 iscode=index(controlcard,'ONE_LETTER')
248 indphi=index(controlcard,'PHI')
249 indback=index(controlcard,'BACK')
250 iranconf=index(controlcard,'RAND_CONF')
251 i2ndstr=index(controlcard,'USE_SEC_PRED')
252 gradout=index(controlcard,'GRADOUT').gt.0
253 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
254 C DISTCHAINMAX become obsolete for periodic boundry condition
255 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
256 C Reading the dimensions of box in x,y,z coordinates
257 call reada(controlcard,'BOXX',boxxsize,100.0d0)
258 call reada(controlcard,'BOXY',boxysize,100.0d0)
259 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
260 c Cutoff range for interactions
261 call reada(controlcard,"R_CUT",r_cut,15.0d0)
262 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
263 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
264 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
265 if (lipthick.gt.0.0d0) then
266 bordliptop=(boxzsize+lipthick)/2.0
267 bordlipbot=bordliptop-lipthick
269 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
270 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
271 buflipbot=bordlipbot+lipbufthick
272 bufliptop=bordliptop-lipbufthick
273 if ((lipbufthick*2.0d0).gt.lipthick)
274 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
276 c write(iout,*) "bordliptop=",bordliptop
277 c write(iout,*) "bordlipbot=",bordlipbot
278 c write(iout,*) "bufliptop=",bufliptop
279 c write(iout,*) "buflipbot=",buflipbot
282 if (me.eq.king .or. .not.out1file )
283 & write (iout,*) "DISTCHAINMAX",distchainmax
285 if(me.eq.king.or..not.out1file)
286 & write (iout,'(2a)') diagmeth(kdiag),
287 & ' routine used to diagonalize matrices.'
290 c--------------------------------------------------------------------------
291 subroutine read_REMDpar
295 implicit real*8 (a-h,o-z)
297 include 'COMMON.IOUNITS'
298 include 'COMMON.TIME1'
301 include 'COMMON.LANGEVIN'
303 include 'COMMON.LANGEVIN.lang0'
305 include 'COMMON.INTERACT'
306 include 'COMMON.NAMES'
308 include 'COMMON.REMD'
309 include 'COMMON.CONTROL'
310 include 'COMMON.SETUP'
312 character*320 controlcard
313 character*3200 controlcard1
314 integer iremd_m_total
316 if(me.eq.king.or..not.out1file)
317 & write (iout,*) "REMD setup"
319 call card_concat(controlcard)
320 call readi(controlcard,"NREP",nrep,3)
321 call readi(controlcard,"NSTEX",nstex,1000)
322 call reada(controlcard,"RETMIN",retmin,10.0d0)
323 call reada(controlcard,"RETMAX",retmax,1000.0d0)
324 mremdsync=(index(controlcard,'SYNC').gt.0)
325 call readi(controlcard,"NSYN",i_sync_step,100)
326 restart1file=(index(controlcard,'REST1FILE').gt.0)
327 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
328 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
329 if(max_cache_traj_use.gt.max_cache_traj)
330 & max_cache_traj_use=max_cache_traj
331 if(me.eq.king.or..not.out1file) then
332 cd if (traj1file) then
333 crc caching is in testing - NTWX is not ignored
334 cd write (iout,*) "NTWX value is ignored"
335 cd write (iout,*) " trajectory is stored to one file by master"
336 cd write (iout,*) " before exchange at NSTEX intervals"
338 write (iout,*) "NREP= ",nrep
339 write (iout,*) "NSTEX= ",nstex
340 write (iout,*) "SYNC= ",mremdsync
341 write (iout,*) "NSYN= ",i_sync_step
342 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
345 if (index(controlcard,'TLIST').gt.0) then
347 call card_concat(controlcard1)
348 read(controlcard1,*) (remd_t(i),i=1,nrep)
349 if(me.eq.king.or..not.out1file)
350 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
353 if (index(controlcard,'MLIST').gt.0) then
355 call card_concat(controlcard1)
356 read(controlcard1,*) (remd_m(i),i=1,nrep)
357 if(me.eq.king.or..not.out1file) then
358 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
361 iremd_m_total=iremd_m_total+remd_m(i)
363 write (iout,*) 'Total number of replicas ',iremd_m_total
366 if(me.eq.king.or..not.out1file)
367 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
370 c--------------------------------------------------------------------------
371 subroutine read_MDpar
375 implicit real*8 (a-h,o-z)
377 include 'COMMON.IOUNITS'
378 include 'COMMON.TIME1'
381 include 'COMMON.LANGEVIN'
383 include 'COMMON.LANGEVIN.lang0'
385 include 'COMMON.INTERACT'
386 include 'COMMON.NAMES'
388 include 'COMMON.SETUP'
389 include 'COMMON.CONTROL'
390 include 'COMMON.SPLITELE'
392 character*320 controlcard
394 call card_concat(controlcard)
395 call readi(controlcard,"NSTEP",n_timestep,1000000)
396 call readi(controlcard,"NTWE",ntwe,100)
397 call readi(controlcard,"NTWX",ntwx,1000)
398 call reada(controlcard,"DT",d_time,1.0d-1)
399 call reada(controlcard,"DVMAX",dvmax,2.0d1)
400 call reada(controlcard,"DAMAX",damax,1.0d1)
401 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
402 call readi(controlcard,"LANG",lang,0)
403 RESPA = index(controlcard,"RESPA") .gt. 0
404 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
405 ntime_split0=ntime_split
406 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
407 ntime_split0=ntime_split
408 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
409 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
410 rest = index(controlcard,"REST").gt.0
411 tbf = index(controlcard,"TBF").gt.0
412 usampl = index(controlcard,"USAMPL").gt.0
414 mdpdb = index(controlcard,"MDPDB").gt.0
415 call reada(controlcard,"T_BATH",t_bath,300.0d0)
416 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
417 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
418 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
419 if (count_reset_moment.eq.0) count_reset_moment=1000000000
420 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
421 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
422 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
423 if (count_reset_vel.eq.0) count_reset_vel=1000000000
424 large = index(controlcard,"LARGE").gt.0
425 print_compon = index(controlcard,"PRINT_COMPON").gt.0
426 rattle = index(controlcard,"RATTLE").gt.0
427 preminim = index(controlcard,"PREMINIM").gt.0
429 dccart=(index(controlcard,'CART').gt.0)
432 c if performing umbrella sampling, fragments constrained are read from the fragment file
438 if(me.eq.king.or..not.out1file) then
440 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
442 write (iout,'(a)') "The units are:"
443 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
444 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
445 & " acceleration: angstrom/(48.9 fs)**2"
446 write (iout,'(a)') "energy: kcal/mol, temperature: K"
448 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
449 write (iout,'(a60,f10.5,a)')
450 & "Initial time step of numerical integration:",d_time,
452 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
454 write (iout,'(2a,i4,a)')
455 & "A-MTS algorithm used; initial time step for fast-varying",
456 & " short-range forces split into",ntime_split," steps."
457 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
458 & r_cut," lambda",rlamb
460 write (iout,'(2a,f10.5)')
461 & "Maximum acceleration threshold to reduce the time step",
462 & "/increase split number:",damax
463 write (iout,'(2a,f10.5)')
464 & "Maximum predicted energy drift to reduce the timestep",
465 & "/increase split number:",edriftmax
466 write (iout,'(a60,f10.5)')
467 & "Maximum velocity threshold to reduce velocities:",dvmax
468 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
469 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
470 if (rattle) write (iout,'(a60)')
471 & "Rattle algorithm used to constrain the virtual bonds"
472 if (preminim .or. iranconf.gt.0) then
474 & "Initial structure will be energy-minimized"
479 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
480 call reada(controlcard,"RWAT",rwat,1.4d0)
481 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
482 surfarea=index(controlcard,"SURFAREA").gt.0
483 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
484 if(me.eq.king.or..not.out1file)then
485 write (iout,'(/a,$)') "Langevin dynamics calculation"
488 & " with direct integration of Langevin equations"
489 else if (lang.eq.2) then
490 write (iout,'(a/)') " with TINKER stochasic MD integrator"
491 else if (lang.eq.3) then
492 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
493 else if (lang.eq.4) then
494 write (iout,'(a/)') " in overdamped mode"
496 write (iout,'(//a,i5)')
497 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
500 write (iout,'(a60,f10.5)') "Temperature:",t_bath
501 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
502 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
503 write (iout,'(a60,f10.5)')
504 & "Scaling factor of the friction forces:",scal_fric
505 if (surfarea) write (iout,'(2a,i10,a)')
506 & "Friction coefficients will be scaled by solvent-accessible",
507 & " surface area every",reset_fricmat," steps."
509 c Calculate friction coefficients and bounds of stochastic forces
510 eta=6*pi*cPoise*etawat
511 if(me.eq.king.or..not.out1file)
512 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
514 gamp=scal_fric*(pstok+rwat)*eta
515 stdfp=dsqrt(2*Rb*t_bath/d_time)
517 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
518 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
520 if(me.eq.king.or..not.out1file)then
521 write (iout,'(/2a/)')
522 & "Radii of site types and friction coefficients and std's of",
523 & " stochastic forces of fully exposed sites"
524 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
526 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
527 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
531 if(me.eq.king.or..not.out1file)then
532 write (iout,'(a)') "Berendsen bath calculation"
533 write (iout,'(a60,f10.5)') "Temperature:",t_bath
534 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
536 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
537 & count_reset_moment," steps"
539 & write (iout,'(a,i10,a)')
540 & "Velocities will be reset at random every",count_reset_vel,
544 if(me.eq.king.or..not.out1file)
545 & write (iout,'(a31)') "Microcanonical mode calculation"
547 if(me.eq.king.or..not.out1file)then
548 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
550 write(iout,*) "MD running with constraints."
551 write(iout,*) "Equilibration time ", eq_time, " mtus."
552 write(iout,*) "Constraining ", nfrag," fragments."
553 write(iout,*) "Length of each fragment, weight and q0:"
555 write (iout,*) "Set of restraints #",iset
557 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
558 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
560 write(iout,*) "constraints between ", npair, "fragments."
561 write(iout,*) "constraint pairs, weights and q0:"
563 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
564 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
566 write(iout,*) "angle constraints within ", nfrag_back,
567 & "backbone fragments."
568 write(iout,*) "fragment, weights:"
570 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
571 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
572 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
575 iset=mod(kolor,nset)+1
578 if(me.eq.king.or..not.out1file)
579 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
582 c------------------------------------------------------------------------------
585 C Read molecular data.
587 implicit real*8 (a-h,o-z)
593 include 'COMMON.IOUNITS'
596 include 'COMMON.INTERACT'
597 include 'COMMON.LOCAL'
598 include 'COMMON.NAMES'
599 include 'COMMON.CHAIN'
600 include 'COMMON.FFIELD'
601 include 'COMMON.SBRIDGE'
602 include 'COMMON.HEADER'
603 include 'COMMON.CONTROL'
604 include 'COMMON.DBASE'
605 include 'COMMON.THREAD'
606 include 'COMMON.CONTACTS'
607 include 'COMMON.TORCNSTR'
608 include 'COMMON.TIME1'
609 include 'COMMON.BOUNDS'
611 include 'COMMON.SETUP'
612 character*4 sequence(maxres)
614 double precision x(maxvar)
615 character*256 pdbfile
616 character*320 weightcard
617 character*80 weightcard_t,ucase
618 dimension itype_pdb(maxres)
619 common /pizda/ itype_pdb
620 logical seq_comp,fail
621 double precision energia(0:n_ene)
627 C Read weights of the subsequent energy terms.
628 call card_concat(weightcard)
629 call reada(weightcard,'WLONG',wlong,1.0D0)
630 call reada(weightcard,'WSC',wsc,wlong)
631 call reada(weightcard,'WSCP',wscp,wlong)
632 call reada(weightcard,'WELEC',welec,1.0D0)
633 call reada(weightcard,'WVDWPP',wvdwpp,welec)
634 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
635 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
636 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
637 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
638 call reada(weightcard,'WTURN3',wturn3,1.0D0)
639 call reada(weightcard,'WTURN4',wturn4,1.0D0)
640 call reada(weightcard,'WTURN6',wturn6,1.0D0)
641 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
642 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
643 call reada(weightcard,'WBOND',wbond,1.0D0)
644 call reada(weightcard,'WTOR',wtor,1.0D0)
645 call reada(weightcard,'WTORD',wtor_d,1.0D0)
646 call reada(weightcard,'WANG',wang,1.0D0)
647 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
648 call reada(weightcard,'SCAL14',scal14,0.4D0)
649 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
650 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
651 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
652 call reada(weightcard,'TEMP0',temp0,300.0d0)
653 call reada(weightcard,'WLT',wliptran,0.0D0)
654 call reada(weightcard,'WSAXS',wsaxs,1.0D0)
655 if (index(weightcard,'SOFT').gt.0) ipot=6
656 C 12/1/95 Added weight for the multi-body term WCORR
657 call reada(weightcard,'WCORRH',wcorr,1.0D0)
658 if (wcorr4.gt.0.0d0) wcorr=wcorr4
679 if(me.eq.king.or..not.out1file)
680 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
681 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
683 10 format (/'Energy-term weights (unscaled):'//
684 & 'WSCC= ',f10.6,' (SC-SC)'/
685 & 'WSCP= ',f10.6,' (SC-p)'/
686 & 'WELEC= ',f10.6,' (p-p electr)'/
687 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
688 & 'WBOND= ',f10.6,' (stretching)'/
689 & 'WANG= ',f10.6,' (bending)'/
690 & 'WSCLOC= ',f10.6,' (SC local)'/
691 & 'WTOR= ',f10.6,' (torsional)'/
692 & 'WTORD= ',f10.6,' (double torsional)'/
693 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
694 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
695 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
696 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
697 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
698 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
699 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
700 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
701 & 'WTURN6= ',f10.6,' (turns, 6th order)')
702 if(me.eq.king.or..not.out1file)then
703 if (wcorr4.gt.0.0d0) then
704 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
705 & 'between contact pairs of peptide groups'
706 write (iout,'(2(a,f5.3/))')
707 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
708 & 'Range of quenching the correlation terms:',2*delt_corr
709 else if (wcorr.gt.0.0d0) then
710 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
711 & 'between contact pairs of peptide groups'
713 write (iout,'(a,f8.3)')
714 & 'Scaling factor of 1,4 SC-p interactions:',scal14
715 write (iout,'(a,f8.3)')
716 & 'General scaling factor of SC-p interactions:',scalscp
718 r0_corr=cutoff_corr-delt_corr
720 aad(i,1)=scalscp*aad(i,1)
721 aad(i,2)=scalscp*aad(i,2)
722 bad(i,1)=scalscp*bad(i,1)
723 bad(i,2)=scalscp*bad(i,2)
725 call rescale_weights(t_bath)
726 if(me.eq.king.or..not.out1file)
727 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
728 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
730 22 format (/'Energy-term weights (scaled):'//
731 & 'WSCC= ',f10.6,' (SC-SC)'/
732 & 'WSCP= ',f10.6,' (SC-p)'/
733 & 'WELEC= ',f10.6,' (p-p electr)'/
734 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
735 & 'WBOND= ',f10.6,' (stretching)'/
736 & 'WANG= ',f10.6,' (bending)'/
737 & 'WSCLOC= ',f10.6,' (SC local)'/
738 & 'WTOR= ',f10.6,' (torsional)'/
739 & 'WTORD= ',f10.6,' (double torsional)'/
740 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
741 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
742 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
743 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
744 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
745 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
746 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
747 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
748 & 'WTURN6= ',f10.6,' (turns, 6th order)')
749 if(me.eq.king.or..not.out1file)
750 & write (iout,*) "Reference temperature for weights calculation:",
752 call reada(weightcard,"D0CM",d0cm,3.78d0)
753 call reada(weightcard,"AKCM",akcm,15.1d0)
754 call reada(weightcard,"AKTH",akth,11.0d0)
755 call reada(weightcard,"AKCT",akct,12.0d0)
756 call reada(weightcard,"V1SS",v1ss,-1.08d0)
757 call reada(weightcard,"V2SS",v2ss,7.61d0)
758 call reada(weightcard,"V3SS",v3ss,13.7d0)
759 call reada(weightcard,"EBR",ebr,-5.50D0)
760 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
762 dyn_ss_mask(i)=.false.
766 dyn_ssbond_ij(i,j)=1.0d300
769 call reada(weightcard,"HT",Ht,0.0D0)
771 ss_depth=ebr/wsc-0.25*eps(1,1)
772 Ht=Ht/wsc-0.25*eps(1,1)
773 akcm=akcm*wstrain/wsc
774 akth=akth*wstrain/wsc
775 akct=akct*wstrain/wsc
776 v1ss=v1ss*wstrain/wsc
777 v2ss=v2ss*wstrain/wsc
778 v3ss=v3ss*wstrain/wsc
780 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
783 if(me.eq.king.or..not.out1file) then
784 write (iout,*) "Parameters of the SS-bond potential:"
785 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
787 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
788 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
789 write (iout,*)" HT",Ht
790 print *,'indpdb=',indpdb,' pdbref=',pdbref
792 if (indpdb.gt.0 .or. pdbref) then
793 read(inp,'(a)') pdbfile
794 if(me.eq.king.or..not.out1file)
795 & write (iout,'(2a)') 'PDB data will be read from file ',
796 & pdbfile(:ilen(pdbfile))
797 open(ipdbin,file=pdbfile,status='old',err=33)
799 33 write (iout,'(a)') 'Error opening PDB file.'
802 c print *,'Begin reading pdb data'
811 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
812 & (crefjlee(j,i+nres),j=1,3)
815 c print *,'Finished reading pdb data'
816 if(me.eq.king.or..not.out1file)
817 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
818 & ' nstart_sup=',nstart_sup
820 itype_pdb(i)=itype(i)
824 nct=nstart_sup+nsup-1
825 call contact(.false.,ncont_ref,icont_ref,co)
828 if(me.eq.king.or..not.out1file)
829 & write(iout,*)'Adding sidechains'
833 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
836 do while (fail.and.nsi.le.maxsi)
837 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
841 c Calculalte only the coordinates of the current sidechain; no need to rebuild
843 call locate_side_chain(i)
844 if(fail) write(iout,*)'Adding sidechain failed for res ',
845 & i,' after ',nsi,' trials'
848 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
852 if (indpdb.eq.0) then
853 C Read sequence if not taken from the pdb file.
855 c print *,'nres=',nres
856 if (iscode.gt.0) then
857 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
859 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
861 C Convert sequence to numeric code
863 itype(i)=rescode(i,sequence(i),iscode)
865 C Assign initial virtual bond lengths
871 vbld(i+nres)=dsc(iabs(itype(i)))
872 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
873 c write (iout,*) "i",i," itype",itype(i),
874 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
878 c print '(20i4)',(itype(i),i=1,nres)
881 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
883 if (itype(i).eq.ntyp1) then
887 else if (iabs(itype(i+1)).ne.20) then
889 else if (iabs(itype(i)).ne.20) then
896 if(me.eq.king.or..not.out1file)then
897 write (iout,*) "ITEL"
899 write (iout,*) i,itype(i),itel(i)
901 print *,'Call Read_Bridge.'
904 C 8/13/98 Set limits to generating the dihedral angles
909 read (inp,*) ndih_constr
910 write (iout,*) "ndish_constr",ndih_constr
911 if (ndih_constr.gt.0) then
913 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
914 if(me.eq.king.or..not.out1file)then
916 & 'There are',ndih_constr,' constraints on phi angles.'
918 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
922 phi0(i)=deg2rad*phi0(i)
923 drange(i)=deg2rad*drange(i)
925 if(me.eq.king.or..not.out1file)
926 & write (iout,*) 'FTORS',ftors
929 phibound(1,ii) = phi0(i)-drange(i)
930 phibound(2,ii) = phi0(i)+drange(i)
937 write (iout,'(a)') 'Boundaries in phi angle sampling:'
939 write (iout,'(a3,i5,2f10.1)')
940 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
946 cd print *,'NNT=',NNT,' NCT=',NCT
947 if (itype(1).eq.ntyp1) nnt=2
948 if (itype(nres).eq.ntyp1) nct=nct-1
950 if(me.eq.king.or..not.out1file)
951 & write (iout,'(a,i3)') 'nsup=',nsup
953 if (nsup.le.(nct-nnt+1)) then
954 do i=0,nct-nnt+1-nsup
955 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
961 & 'Error - sequences to be superposed do not match.'
964 do i=0,nsup-(nct-nnt+1)
965 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
967 nstart_sup=nstart_sup+i
973 & 'Error - sequences to be superposed do not match.'
976 if (nsup.eq.0) nsup=nct-nnt
977 if (nstart_sup.eq.0) nstart_sup=nnt
978 if (nstart_seq.eq.0) nstart_seq=nnt
979 if(me.eq.king.or..not.out1file)
980 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
981 & ' nstart_seq=',nstart_seq
983 c--- Zscore rms -------
984 if (nz_start.eq.0) nz_start=nnt
985 if (nz_end.eq.0 .and. nsup.gt.0) then
987 else if (nz_end.eq.0) then
990 if(me.eq.king.or..not.out1file)then
991 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
992 write (iout,*) 'IZ_SC=',iz_sc
994 c----------------------
997 if (.not.pdbref) then
998 call read_angles(inp,*38)
1000 38 write (iout,'(a)') 'Error reading reference structure.'
1002 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1003 stop 'Error reading reference structure'
1005 39 call chainbuild_extconf
1007 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1014 cref(j,i,kkk)=c(j,i)
1017 call contact(.true.,ncont_ref,icont_ref,co)
1021 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1023 if (constr_dist.gt.0) call read_dist_constr
1024 write (iout,*) "After read_dist_constr nhpb",nhpb
1025 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1026 if(me.eq.king.or..not.out1file)
1027 & write (iout,*) 'Contact order:',co
1029 if(me.eq.king.or..not.out1file)
1030 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1033 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1035 if(me.eq.king.or..not.out1file)
1036 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1037 & icont_ref(1,i),' ',
1038 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1041 write (iout,*) "calling read_saxs_consrtr",nsaxs
1042 if (nsaxs.gt.0) call read_saxs_constr
1044 if (constr_homology.gt.0) then
1045 call read_constr_homology
1046 if (indpdb.gt.0 .or. pdbref) then
1049 c(j,i)=crefjlee(j,i)
1050 cref(j,i,1)=crefjlee(j,i)
1055 write (iout,*) "Array C"
1057 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1058 & (c(j,i+nres),j=1,3)
1060 write (iout,*) "Array Cref"
1062 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),
1063 & (cref(j,i+nres,1),j=1,3)
1066 call int_from_cart1(.false.)
1067 call sc_loc_geom(.false.)
1069 thetaref(i)=theta(i)
1074 dc(j,i)=c(j,i+1)-c(j,i)
1075 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1080 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1081 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1089 if (nhpb.gt.0) call hpb_partition
1090 if (peak.gt.0) call NMRpeak_partition
1091 c write (iout,*) "After read_dist_constr nhpb",nhpb
1093 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1094 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1095 & modecalc.ne.10) then
1096 C If input structure hasn't been supplied from the PDB file read or generate
1098 if (iranconf.eq.0 .and. .not. extconf) then
1099 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1100 & write (iout,'(a)') 'Initial geometry will be read in.'
1102 read(inp,'(8f10.5)',end=36,err=36)
1103 & ((c(l,k),l=1,3),k=1,nres),
1104 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1105 write (iout,*) "Exit READ_CART"
1106 write (iout,'(8f10.5)')
1107 & ((c(l,k),l=1,3),k=1,nres),
1108 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1109 call int_from_cart1(.true.)
1110 write (iout,*) "Finish INT_TO_CART"
1113 dc(j,i)=c(j,i+1)-c(j,i)
1114 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1118 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1120 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1121 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1127 write (iout,*) "Calling read_ang"
1128 call read_angles(inp,*36)
1129 write (iout,*) "Calling chainbuild"
1130 call chainbuild_extconf
1133 36 write (iout,'(a)') 'Error reading angle file.'
1135 call mpi_finalize( MPI_COMM_WORLD,IERR )
1137 stop 'Error reading angle file.'
1139 else if (extconf) then
1140 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1141 & write (iout,'(a)') 'Extended chain initial geometry.'
1143 theta(i)=90d0*deg2rad
1146 phi(i)=180d0*deg2rad
1149 alph(i)=110d0*deg2rad
1152 omeg(i)=-120d0*deg2rad
1153 if (itype(i).le.0) omeg(i)=-omeg(i)
1155 c from old chainbuild
1157 C Define the origin and orientation of the coordinate system and locate the
1158 C first three CA's and SC(2).
1162 * Build the alpha-carbon chain.
1165 call locate_next_res(i)
1168 C First and last SC must coincide with the corresponding CA.
1172 dc_norm(j,nres+1)=0.0D0
1173 dc(j,nres+nres)=0.0D0
1174 dc_norm(j,nres+nres)=0.0D0
1176 c(j,nres+nres)=c(j,nres)
1179 C Define the origin and orientation of the coordinate system and locate the
1180 C first three CA's and SC(2).
1184 * Build the alpha-carbon chain.
1187 call locate_next_res(i)
1190 C First and last SC must coincide with the corresponding CA.
1194 dc_norm(j,nres+1)=0.0D0
1195 dc(j,nres+nres)=0.0D0
1196 dc_norm(j,nres+nres)=0.0D0
1198 c(j,nres+nres)=c(j,nres)
1203 if(me.eq.king.or..not.out1file)
1204 & write (iout,'(a)') 'Random-generated initial geometry.'
1208 if (me.eq.king .or. fg_rank.eq.0 .and. (
1209 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1213 call gen_rand_conf(itmp,*30)
1215 30 write (iout,*) 'Failed to generate random conformation',
1216 & ', itrial=',itrial
1217 write (*,*) 'Processor:',me,
1218 & ' Failed to generate random conformation',
1228 write (iout,'(a,i3,a)') 'Processor:',me,
1229 & ' error in generating random conformation.'
1230 write (*,'(a,i3,a)') 'Processor:',me,
1231 & ' error in generating random conformation.'
1234 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1239 & ' error in generating random conformation.'
1244 elseif (modecalc.eq.4) then
1245 read (inp,'(a)') intinname
1246 open (intin,file=intinname,status='old',err=333)
1247 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1248 & write (iout,'(a)') 'intinname',intinname
1249 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1251 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1253 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1255 stop 'Error opening angle file.'
1259 C Generate distance constraints, if the PDB structure is to be regularized.
1260 if (nthread.gt.0) then
1261 call read_threadbase
1264 if (me.eq.king .or. .not. out1file)
1266 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1267 write (iout,'(/a,i3,a)')
1268 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1269 write (iout,'(20i4)') (iss(i),i=1,ns)
1271 write(iout,*)"Running with dynamic disulfide-bond formation"
1273 write (iout,'(/a/)') 'Pre-formed links are:'
1279 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1280 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1286 if (ns.gt.0.and.dyn_ss) then
1290 forcon(i-nss)=forcon(i)
1297 dyn_ss_mask(iss(i))=.true.
1300 if (i2ndstr.gt.0) call secstrp2dihc
1301 c call geom_to_var(nvar,x)
1302 c call etotal(energia(0))
1303 c call enerprint(energia(0))
1304 c call briefout(0,etot)
1306 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1307 cd write (iout,'(a)') 'Variable list:'
1308 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1310 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1311 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1312 & 'Processor',myrank,': end reading molecular data.'
1316 c--------------------------------------------------------------------------
1317 logical function seq_comp(itypea,itypeb,length)
1319 integer length,itypea(length),itypeb(length)
1322 if (itypea(i).ne.itypeb(i)) then
1330 c-----------------------------------------------------------------------------
1331 subroutine read_bridge
1332 C Read information about disulfide bridges.
1333 implicit real*8 (a-h,o-z)
1334 include 'DIMENSIONS'
1338 include 'COMMON.IOUNITS'
1339 include 'COMMON.GEO'
1340 include 'COMMON.VAR'
1341 include 'COMMON.INTERACT'
1342 include 'COMMON.LOCAL'
1343 include 'COMMON.NAMES'
1344 include 'COMMON.CHAIN'
1345 include 'COMMON.FFIELD'
1346 include 'COMMON.SBRIDGE'
1347 include 'COMMON.HEADER'
1348 include 'COMMON.CONTROL'
1349 include 'COMMON.DBASE'
1350 include 'COMMON.THREAD'
1351 include 'COMMON.TIME1'
1352 include 'COMMON.SETUP'
1353 C Read bridging residues.
1354 read (inp,*) ns,(iss(i),i=1,ns)
1356 if(me.eq.king.or..not.out1file)
1357 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1358 C Check whether the specified bridging residues are cystines.
1360 if (itype(iss(i)).ne.1) then
1361 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1362 & 'Do you REALLY think that the residue ',
1363 & restyp(itype(iss(i))),i,
1364 & ' can form a disulfide bridge?!!!'
1365 write (*,'(2a,i3,a)')
1366 & 'Do you REALLY think that the residue ',
1367 & restyp(itype(iss(i))),i,
1368 & ' can form a disulfide bridge?!!!'
1370 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1375 C Read preformed bridges.
1377 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1379 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1382 C Check if the residues involved in bridges are in the specified list of
1383 C bridging residues.
1386 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1387 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1388 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1389 & ' contains residues present in other pairs.'
1390 write (*,'(a,i3,a)') 'Disulfide pair',i,
1391 & ' contains residues present in other pairs.'
1393 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1399 if (ihpb(i).eq.iss(j)) goto 10
1401 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1404 if (jhpb(i).eq.iss(j)) goto 20
1406 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1412 ihpb(i)=ihpb(i)+nres
1413 jhpb(i)=jhpb(i)+nres
1419 c----------------------------------------------------------------------------
1420 subroutine read_x(kanal,*)
1421 implicit real*8 (a-h,o-z)
1422 include 'DIMENSIONS'
1423 include 'COMMON.GEO'
1424 include 'COMMON.VAR'
1425 include 'COMMON.CHAIN'
1426 include 'COMMON.IOUNITS'
1427 include 'COMMON.CONTROL'
1428 include 'COMMON.LOCAL'
1429 include 'COMMON.INTERACT'
1430 c Read coordinates from input
1432 read(kanal,'(8f10.5)',end=10,err=10)
1433 & ((c(l,k),l=1,3),k=1,nres),
1434 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1437 c(j,2*nres)=c(j,nres)
1439 call int_from_cart1(.false.)
1442 dc(j,i)=c(j,i+1)-c(j,i)
1443 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1447 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1449 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1450 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1458 c----------------------------------------------------------------------------
1459 subroutine read_threadbase
1460 implicit real*8 (a-h,o-z)
1461 include 'DIMENSIONS'
1462 include 'COMMON.IOUNITS'
1463 include 'COMMON.GEO'
1464 include 'COMMON.VAR'
1465 include 'COMMON.INTERACT'
1466 include 'COMMON.LOCAL'
1467 include 'COMMON.NAMES'
1468 include 'COMMON.CHAIN'
1469 include 'COMMON.FFIELD'
1470 include 'COMMON.SBRIDGE'
1471 include 'COMMON.HEADER'
1472 include 'COMMON.CONTROL'
1473 include 'COMMON.DBASE'
1474 include 'COMMON.THREAD'
1475 include 'COMMON.TIME1'
1476 C Read pattern database for threading.
1477 read (icbase,*) nseq
1479 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1480 & nres_base(2,i),nres_base(3,i)
1481 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1483 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1484 c & nres_base(2,i),nres_base(3,i)
1485 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1489 if (weidis.eq.0.0D0) weidis=0.1D0
1498 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1499 write (iout,'(a,i5)') 'nexcl: ',nexcl
1500 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1503 c------------------------------------------------------------------------------
1504 subroutine setup_var
1505 implicit real*8 (a-h,o-z)
1506 include 'DIMENSIONS'
1507 include 'COMMON.IOUNITS'
1508 include 'COMMON.GEO'
1509 include 'COMMON.VAR'
1510 include 'COMMON.INTERACT'
1511 include 'COMMON.LOCAL'
1512 include 'COMMON.NAMES'
1513 include 'COMMON.CHAIN'
1514 include 'COMMON.FFIELD'
1515 include 'COMMON.SBRIDGE'
1516 include 'COMMON.HEADER'
1517 include 'COMMON.CONTROL'
1518 include 'COMMON.DBASE'
1519 include 'COMMON.THREAD'
1520 include 'COMMON.TIME1'
1521 C Set up variable list.
1527 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1529 ialph(i,1)=nvar+nside
1533 if (indphi.gt.0) then
1535 else if (indback.gt.0) then
1540 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1543 c----------------------------------------------------------------------------
1544 subroutine gen_dist_constr
1545 C Generate CA distance constraints.
1546 implicit real*8 (a-h,o-z)
1547 include 'DIMENSIONS'
1548 include 'COMMON.IOUNITS'
1549 include 'COMMON.GEO'
1550 include 'COMMON.VAR'
1551 include 'COMMON.INTERACT'
1552 include 'COMMON.LOCAL'
1553 include 'COMMON.NAMES'
1554 include 'COMMON.CHAIN'
1555 include 'COMMON.FFIELD'
1556 include 'COMMON.SBRIDGE'
1557 include 'COMMON.HEADER'
1558 include 'COMMON.CONTROL'
1559 include 'COMMON.DBASE'
1560 include 'COMMON.THREAD'
1561 include 'COMMON.TIME1'
1562 dimension itype_pdb(maxres)
1563 common /pizda/ itype_pdb
1565 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1566 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1567 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1569 do i=nstart_sup,nstart_sup+nsup-1
1570 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1571 cd & ' seq_pdb', restyp(itype_pdb(i))
1572 do j=i+2,nstart_sup+nsup-1
1574 ihpb(nhpb)=i+nstart_seq-nstart_sup
1575 jhpb(nhpb)=j+nstart_seq-nstart_sup
1577 dhpb(nhpb)=dist(i,j)
1580 cd write (iout,'(a)') 'Distance constraints:'
1585 cd if (ii.gt.nres) then
1590 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1591 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1592 cd & dhpb(i),forcon(i)
1596 c----------------------------------------------------------------------------
1598 implicit real*8 (a-h,o-z)
1599 include 'DIMENSIONS'
1600 include 'COMMON.MAP'
1601 include 'COMMON.IOUNITS'
1602 character*3 angid(4) /'THE','PHI','ALP','OME'/
1603 character*80 mapcard,ucase
1605 read (inp,'(a)') mapcard
1606 mapcard=ucase(mapcard)
1607 if (index(mapcard,'PHI').gt.0) then
1609 else if (index(mapcard,'THE').gt.0) then
1611 else if (index(mapcard,'ALP').gt.0) then
1613 else if (index(mapcard,'OME').gt.0) then
1616 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1617 stop 'Error - illegal variable spec in MAP card.'
1619 call readi (mapcard,'RES1',res1(imap),0)
1620 call readi (mapcard,'RES2',res2(imap),0)
1621 if (res1(imap).eq.0) then
1622 res1(imap)=res2(imap)
1623 else if (res2(imap).eq.0) then
1624 res2(imap)=res1(imap)
1626 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1628 & 'Error - illegal definition of variable group in MAP.'
1629 stop 'Error - illegal definition of variable group in MAP.'
1631 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1632 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1633 call readi(mapcard,'NSTEP',nstep(imap),0)
1634 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1636 & 'Illegal boundary and/or step size specification in MAP.'
1637 stop 'Illegal boundary and/or step size specification in MAP.'
1642 c----------------------------------------------------------------------------
1644 implicit real*8 (a-h,o-z)
1645 include 'DIMENSIONS'
1646 include 'COMMON.IOUNITS'
1647 include 'COMMON.GEO'
1648 include 'COMMON.CSA'
1649 include 'COMMON.BANK'
1650 include 'COMMON.CONTROL'
1652 character*620 mcmcard
1653 call card_concat(mcmcard)
1655 call readi(mcmcard,'NCONF',nconf,50)
1656 call readi(mcmcard,'NADD',nadd,0)
1657 call readi(mcmcard,'JSTART',jstart,1)
1658 call readi(mcmcard,'JEND',jend,1)
1659 call readi(mcmcard,'NSTMAX',nstmax,500000)
1660 call readi(mcmcard,'N0',n0,1)
1661 call readi(mcmcard,'N1',n1,6)
1662 call readi(mcmcard,'N2',n2,4)
1663 call readi(mcmcard,'N3',n3,0)
1664 call readi(mcmcard,'N4',n4,0)
1665 call readi(mcmcard,'N5',n5,0)
1666 call readi(mcmcard,'N6',n6,10)
1667 call readi(mcmcard,'N7',n7,0)
1668 call readi(mcmcard,'N8',n8,0)
1669 call readi(mcmcard,'N9',n9,0)
1670 call readi(mcmcard,'N14',n14,0)
1671 call readi(mcmcard,'N15',n15,0)
1672 call readi(mcmcard,'N16',n16,0)
1673 call readi(mcmcard,'N17',n17,0)
1674 call readi(mcmcard,'N18',n18,0)
1676 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1678 call readi(mcmcard,'NDIFF',ndiff,2)
1679 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1680 call readi(mcmcard,'IS1',is1,1)
1681 call readi(mcmcard,'IS2',is2,8)
1682 call readi(mcmcard,'NRAN0',nran0,4)
1683 call readi(mcmcard,'NRAN1',nran1,2)
1684 call readi(mcmcard,'IRR',irr,1)
1685 call readi(mcmcard,'NSEED',nseed,20)
1686 call readi(mcmcard,'NTOTAL',ntotal,10000)
1687 call reada(mcmcard,'CUT1',cut1,2.0d0)
1688 call reada(mcmcard,'CUT2',cut2,5.0d0)
1689 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1690 call readi(mcmcard,'ICMAX',icmax,3)
1691 call readi(mcmcard,'IRESTART',irestart,0)
1692 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1695 call reada(mcmcard,'DELE',dele,20.0d0)
1696 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1697 call readi(mcmcard,'IREF',iref,0)
1698 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1699 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1700 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1701 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1702 write (iout,*) "NCONF_IN",nconf_in
1705 c----------------------------------------------------------------------------
1706 cfmc subroutine mcmfread
1707 cfmc implicit real*8 (a-h,o-z)
1708 cfmc include 'DIMENSIONS'
1709 cfmc include 'COMMON.MCMF'
1710 cfmc include 'COMMON.IOUNITS'
1711 cfmc include 'COMMON.GEO'
1712 cfmc character*80 ucase
1713 cfmc character*620 mcmcard
1714 cfmc call card_concat(mcmcard)
1716 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1717 cfmc write(iout,*)'MAXRANT=',maxrant
1718 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1719 cfmc write(iout,*)'MAXFAM=',maxfam
1720 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1721 cfmc write(iout,*)'NNET1=',nnet1
1722 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1723 cfmc write(iout,*)'NNET2=',nnet2
1724 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1725 cfmc write(iout,*)'NNET3=',nnet3
1726 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1727 cfmc write(iout,*)'ILASTT=',ilastt
1728 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1729 cfmc write(iout,*)'MAXSTR=',maxstr
1730 cfmc maxstr_f=maxstr/maxfam
1731 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1732 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1733 cfmc write(iout,*)'NMCMF=',nmcmf
1734 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1735 cfmc write(iout,*)'IFOCUS=',ifocus
1736 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1737 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1738 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1739 cfmc write(iout,*)'INTPRT=',intprt
1740 cfmc call readi(mcmcard,'IPRT',iprt,100)
1741 cfmc write(iout,*)'IPRT=',iprt
1742 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1743 cfmc write(iout,*)'IMAXTR=',imaxtr
1744 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1745 cfmc write(iout,*)'MAXEVEN=',maxeven
1746 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1747 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1748 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1749 cfmc write(iout,*)'INIMIN=',inimin
1750 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1751 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1752 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1753 cfmc write(iout,*)'NTHREAD=',nthread
1754 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1755 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1756 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1757 cfmc write(iout,*)'MAXPERT=',maxpert
1758 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1759 cfmc write(iout,*)'IRMSD=',irmsd
1760 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1761 cfmc write(iout,*)'DENEMIN=',denemin
1762 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1763 cfmc write(iout,*)'RCUT1S=',rcut1s
1764 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1765 cfmc write(iout,*)'RCUT1E=',rcut1e
1766 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1767 cfmc write(iout,*)'RCUT2S=',rcut2s
1768 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1769 cfmc write(iout,*)'RCUT2E=',rcut2e
1770 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1771 cfmc write(iout,*)'DPERT1=',d_pert1
1772 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1773 cfmc write(iout,*)'DPERT1A=',d_pert1a
1774 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1775 cfmc write(iout,*)'DPERT2=',d_pert2
1776 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1777 cfmc write(iout,*)'DPERT2A=',d_pert2a
1778 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1779 cfmc write(iout,*)'DPERT2B=',d_pert2b
1780 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1781 cfmc write(iout,*)'DPERT2C=',d_pert2c
1782 cfmc d_pert1=deg2rad*d_pert1
1783 cfmc d_pert1a=deg2rad*d_pert1a
1784 cfmc d_pert2=deg2rad*d_pert2
1785 cfmc d_pert2a=deg2rad*d_pert2a
1786 cfmc d_pert2b=deg2rad*d_pert2b
1787 cfmc d_pert2c=deg2rad*d_pert2c
1788 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1789 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1790 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1791 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1792 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1793 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1794 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1795 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1796 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1797 cfmc write(iout,*)'RCUTINI=',rcutini
1798 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1799 cfmc write(iout,*)'GRAT=',grat
1800 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1801 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1805 c----------------------------------------------------------------------------
1807 implicit real*8 (a-h,o-z)
1808 include 'DIMENSIONS'
1809 include 'COMMON.MCM'
1810 include 'COMMON.MCE'
1811 include 'COMMON.IOUNITS'
1813 character*320 mcmcard
1814 call card_concat(mcmcard)
1815 call readi(mcmcard,'MAXACC',maxacc,100)
1816 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1817 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1818 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1819 call readi(mcmcard,'MAXREPM',maxrepm,200)
1820 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1821 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1822 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1823 call reada(mcmcard,'E_UP',e_up,5.0D0)
1824 call reada(mcmcard,'DELTE',delte,0.1D0)
1825 call readi(mcmcard,'NSWEEP',nsweep,5)
1826 call readi(mcmcard,'NSTEPH',nsteph,0)
1827 call readi(mcmcard,'NSTEPC',nstepc,0)
1828 call reada(mcmcard,'TMIN',tmin,298.0D0)
1829 call reada(mcmcard,'TMAX',tmax,298.0D0)
1830 call readi(mcmcard,'NWINDOW',nwindow,0)
1831 call readi(mcmcard,'PRINT_MC',print_mc,0)
1832 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1833 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1834 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1835 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1836 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1837 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1838 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1839 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1840 if (nwindow.gt.0) then
1841 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1843 winlen(i)=winend(i)-winstart(i)+1
1846 if (tmax.lt.tmin) tmax=tmin
1847 if (tmax.eq.tmin) then
1851 if (nstepc.gt.0 .and. nsteph.gt.0) then
1852 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1853 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1855 C Probabilities of different move types
1856 sumpro_type(0)=0.0D0
1857 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1858 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1859 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1860 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1861 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1862 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1863 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1865 print *,'i',i,' sumprotype',sumpro_type(i)
1866 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1867 print *,'i',i,' sumprotype',sumpro_type(i)
1871 c----------------------------------------------------------------------------
1872 subroutine read_minim
1873 implicit real*8 (a-h,o-z)
1874 include 'DIMENSIONS'
1875 include 'COMMON.MINIM'
1876 include 'COMMON.IOUNITS'
1877 include 'COMMON.CONTROL'
1878 include 'COMMON.SETUP'
1880 character*320 minimcard
1881 call card_concat(minimcard)
1882 call readi(minimcard,'MAXMIN',maxmin,2000)
1883 call readi(minimcard,'MAXFUN',maxfun,5000)
1884 call readi(minimcard,'MINMIN',minmin,maxmin)
1885 call readi(minimcard,'MINFUN',minfun,maxmin)
1886 call reada(minimcard,'TOLF',tolf,1.0D-2)
1887 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1888 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1889 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1890 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1892 if (.not. out1file .or. me.eq.king) then
1894 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1895 & 'Options in energy minimization:'
1896 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1897 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1898 & 'MinMin:',MinMin,' MinFun:',MinFun,
1899 & ' TolF:',TolF,' RTolF:',RTolF
1905 c----------------------------------------------------------------------------
1906 subroutine read_angles(kanal,*)
1907 implicit real*8 (a-h,o-z)
1908 include 'DIMENSIONS'
1909 include 'COMMON.GEO'
1910 include 'COMMON.VAR'
1911 include 'COMMON.CHAIN'
1912 include 'COMMON.IOUNITS'
1913 include 'COMMON.CONTROL'
1914 c Read angles from input
1916 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1917 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1918 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1919 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1922 c 9/7/01 avoid 180 deg valence angle
1923 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1925 theta(i)=deg2rad*theta(i)
1926 phi(i)=deg2rad*phi(i)
1927 alph(i)=deg2rad*alph(i)
1928 omeg(i)=deg2rad*omeg(i)
1933 c----------------------------------------------------------------------------
1934 subroutine reada(rekord,lancuch,wartosc,default)
1936 character*(*) rekord,lancuch
1937 double precision wartosc,default
1940 iread=index(rekord,lancuch)
1941 if (iread.eq.0) then
1945 iread=iread+ilen(lancuch)+1
1946 read (rekord(iread:),*,err=10,end=10) wartosc
1951 c----------------------------------------------------------------------------
1952 subroutine readi(rekord,lancuch,wartosc,default)
1954 character*(*) rekord,lancuch
1955 integer wartosc,default
1958 iread=index(rekord,lancuch)
1959 if (iread.eq.0) then
1963 iread=iread+ilen(lancuch)+1
1964 read (rekord(iread:),*,err=10,end=10) wartosc
1969 c----------------------------------------------------------------------------
1970 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1973 integer tablica(dim),default
1974 character*(*) rekord,lancuch
1981 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1982 if (iread.eq.0) return
1983 iread=iread+ilen(lancuch)+1
1984 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1987 c----------------------------------------------------------------------------
1988 subroutine multreada(rekord,lancuch,tablica,dim,default)
1991 double precision tablica(dim),default
1992 character*(*) rekord,lancuch
1999 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2000 if (iread.eq.0) return
2001 iread=iread+ilen(lancuch)+1
2002 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2005 c----------------------------------------------------------------------------
2006 subroutine openunits
2007 implicit real*8 (a-h,o-z)
2008 include 'DIMENSIONS'
2011 character*16 form,nodename
2014 include 'COMMON.SETUP'
2015 include 'COMMON.IOUNITS'
2017 include 'COMMON.CONTROL'
2018 integer lenpre,lenpot,ilen,lentmp
2020 character*3 out1file_text,ucase
2023 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2024 call getenv_loc("PREFIX",prefix)
2026 call getenv_loc("POT",pot)
2027 call getenv_loc("DIRTMP",tmpdir)
2028 call getenv_loc("CURDIR",curdir)
2029 call getenv_loc("OUT1FILE",out1file_text)
2030 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2031 out1file_text=ucase(out1file_text)
2032 if (out1file_text(1:1).eq."Y") then
2035 out1file=fg_rank.gt.0
2040 if (lentmp.gt.0) then
2041 write (*,'(80(1h!))')
2042 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2043 write (*,'(80(1h!))')
2044 write (*,*)"All output files will be on node /tmp directory."
2046 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2047 if (me.eq.king) then
2048 write (*,*) "The master node is ",nodename
2049 else if (fg_rank.eq.0) then
2050 write (*,*) "I am the CG slave node ",nodename
2052 write (*,*) "I am the FG slave node ",nodename
2055 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2056 lenpre = lentmp+lenpre+1
2058 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2059 C Get the names and open the input files
2060 #if defined(WINIFL) || defined(WINPGI)
2061 open(1,file=pref_orig(:ilen(pref_orig))//
2062 & '.inp',status='old',readonly,shared)
2063 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2064 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2065 C Get parameter filenames and open the parameter files.
2066 call getenv_loc('BONDPAR',bondname)
2067 open (ibond,file=bondname,status='old',readonly,shared)
2068 call getenv_loc('THETPAR',thetname)
2069 open (ithep,file=thetname,status='old',readonly,shared)
2071 call getenv_loc('THETPARPDB',thetname_pdb)
2072 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2074 call getenv_loc('ROTPAR',rotname)
2075 open (irotam,file=rotname,status='old',readonly,shared)
2077 call getenv_loc('ROTPARPDB',rotname_pdb)
2078 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2080 call getenv_loc('TORPAR',torname)
2081 open (itorp,file=torname,status='old',readonly,shared)
2082 call getenv_loc('TORDPAR',tordname)
2083 open (itordp,file=tordname,status='old',readonly,shared)
2084 call getenv_loc('FOURIER',fouriername)
2085 open (ifourier,file=fouriername,status='old',readonly,shared)
2086 call getenv_loc('ELEPAR',elename)
2087 open (ielep,file=elename,status='old',readonly,shared)
2088 call getenv_loc('SIDEPAR',sidename)
2089 open (isidep,file=sidename,status='old',readonly,shared)
2090 #elif (defined CRAY) || (defined AIX)
2091 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2093 c print *,"Processor",myrank," opened file 1"
2094 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2095 c print *,"Processor",myrank," opened file 9"
2096 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2097 C Get parameter filenames and open the parameter files.
2098 call getenv_loc('BONDPAR',bondname)
2099 open (ibond,file=bondname,status='old',action='read')
2100 c print *,"Processor",myrank," opened file IBOND"
2101 call getenv_loc('THETPAR',thetname)
2102 open (ithep,file=thetname,status='old',action='read')
2103 c print *,"Processor",myrank," opened file ITHEP"
2105 call getenv_loc('THETPARPDB',thetname_pdb)
2106 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2108 call getenv_loc('ROTPAR',rotname)
2109 open (irotam,file=rotname,status='old',action='read')
2110 c print *,"Processor",myrank," opened file IROTAM"
2112 call getenv_loc('ROTPARPDB',rotname_pdb)
2113 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2115 call getenv_loc('TORPAR',torname)
2116 open (itorp,file=torname,status='old',action='read')
2117 c print *,"Processor",myrank," opened file ITORP"
2118 call getenv_loc('TORDPAR',tordname)
2119 open (itordp,file=tordname,status='old',action='read')
2120 c print *,"Processor",myrank," opened file ITORDP"
2121 call getenv_loc('SCCORPAR',sccorname)
2122 open (isccor,file=sccorname,status='old',action='read')
2123 c print *,"Processor",myrank," opened file ISCCOR"
2124 call getenv_loc('FOURIER',fouriername)
2125 open (ifourier,file=fouriername,status='old',action='read')
2126 c print *,"Processor",myrank," opened file IFOURIER"
2127 call getenv_loc('ELEPAR',elename)
2128 open (ielep,file=elename,status='old',action='read')
2129 c print *,"Processor",myrank," opened file IELEP"
2130 call getenv_loc('SIDEPAR',sidename)
2131 open (isidep,file=sidename,status='old',action='read')
2132 c print *,"Processor",myrank," opened file ISIDEP"
2133 c print *,"Processor",myrank," opened parameter files"
2135 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2136 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2137 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2138 C Get parameter filenames and open the parameter files.
2139 call getenv_loc('BONDPAR',bondname)
2140 open (ibond,file=bondname,status='old')
2141 call getenv_loc('THETPAR',thetname)
2142 open (ithep,file=thetname,status='old')
2144 call getenv_loc('THETPARPDB',thetname_pdb)
2145 open (ithep_pdb,file=thetname_pdb,status='old')
2147 call getenv_loc('ROTPAR',rotname)
2148 open (irotam,file=rotname,status='old')
2150 call getenv_loc('ROTPARPDB',rotname_pdb)
2151 open (irotam_pdb,file=rotname_pdb,status='old')
2153 call getenv_loc('TORPAR',torname)
2154 open (itorp,file=torname,status='old')
2155 call getenv_loc('TORDPAR',tordname)
2156 open (itordp,file=tordname,status='old')
2157 call getenv_loc('SCCORPAR',sccorname)
2158 open (isccor,file=sccorname,status='old')
2159 call getenv_loc('FOURIER',fouriername)
2160 open (ifourier,file=fouriername,status='old')
2161 call getenv_loc('ELEPAR',elename)
2162 open (ielep,file=elename,status='old')
2163 call getenv_loc('SIDEPAR',sidename)
2164 open (isidep,file=sidename,status='old')
2165 call getenv_loc('LIPTRANPAR',liptranname)
2166 open (iliptranpar,file=liptranname,status='old')
2168 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2170 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2171 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2172 C Get parameter filenames and open the parameter files.
2173 call getenv_loc('BONDPAR',bondname)
2174 open (ibond,file=bondname,status='old',readonly)
2175 call getenv_loc('THETPAR',thetname)
2176 open (ithep,file=thetname,status='old',readonly)
2177 call getenv_loc('ROTPAR',rotname)
2178 open (irotam,file=rotname,status='old',readonly)
2179 call getenv_loc('TORPAR',torname)
2180 open (itorp,file=torname,status='old',readonly)
2181 call getenv_loc('TORDPAR',tordname)
2182 open (itordp,file=tordname,status='old',readonly)
2183 call getenv_loc('SCCORPAR',sccorname)
2184 open (isccor,file=sccorname,status='old',readonly)
2186 call getenv_loc('THETPARPDB',thetname_pdb)
2187 c print *,"thetname_pdb ",thetname_pdb
2188 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2189 c print *,ithep_pdb," opened"
2191 call getenv_loc('FOURIER',fouriername)
2192 open (ifourier,file=fouriername,status='old',readonly)
2193 call getenv_loc('ELEPAR',elename)
2194 open (ielep,file=elename,status='old',readonly)
2195 call getenv_loc('SIDEPAR',sidename)
2196 open (isidep,file=sidename,status='old',readonly)
2197 call getenv_loc('LIPTRANPAR',liptranname)
2198 open (iliptranpar,file=liptranname,status='old',readonly)
2200 call getenv_loc('ROTPARPDB',rotname_pdb)
2201 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2206 C 8/9/01 In the newest version SCp interaction constants are read from a file
2207 C Use -DOLDSCP to use hard-coded constants instead.
2209 call getenv_loc('SCPPAR',scpname)
2210 #if defined(WINIFL) || defined(WINPGI)
2211 open (iscpp,file=scpname,status='old',readonly,shared)
2212 #elif (defined CRAY) || (defined AIX)
2213 open (iscpp,file=scpname,status='old',action='read')
2215 open (iscpp,file=scpname,status='old')
2217 open (iscpp,file=scpname,status='old',readonly)
2220 call getenv_loc('PATTERN',patname)
2221 #if defined(WINIFL) || defined(WINPGI)
2222 open (icbase,file=patname,status='old',readonly,shared)
2223 #elif (defined CRAY) || (defined AIX)
2224 open (icbase,file=patname,status='old',action='read')
2226 open (icbase,file=patname,status='old')
2228 open (icbase,file=patname,status='old',readonly)
2231 C Open output file only for CG processes
2232 c print *,"Processor",myrank," fg_rank",fg_rank
2233 if (fg_rank.eq.0) then
2235 if (nodes.eq.1) then
2238 npos = dlog10(dfloat(nodes-1))+1
2240 if (npos.lt.3) npos=3
2241 write (liczba,'(i1)') npos
2242 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2244 write (liczba,form) me
2245 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2246 & liczba(:ilen(liczba))
2247 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2249 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2251 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2252 & liczba(:ilen(liczba))//'.mol2'
2253 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2254 & liczba(:ilen(liczba))//'.stat'
2256 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2257 & //liczba(:ilen(liczba))//'.stat')
2258 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2261 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2262 & liczba(:ilen(liczba))//'.const'
2267 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2268 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2269 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2270 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2271 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2273 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2275 rest2name=prefix(:ilen(prefix))//'.rst'
2277 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2280 #if defined(AIX) || defined(PGI)
2281 if (me.eq.king .or. .not. out1file)
2282 & open(iout,file=outname,status='unknown')
2284 if (fg_rank.gt.0) then
2285 write (liczba,'(i3.3)') myrank/nfgtasks
2286 write (ll,'(bz,i3.3)') fg_rank
2287 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2292 open(igeom,file=intname,status='unknown',position='append')
2293 open(ipdb,file=pdbname,status='unknown')
2294 open(imol2,file=mol2name,status='unknown')
2295 open(istat,file=statname,status='unknown',position='append')
2297 c1out open(iout,file=outname,status='unknown')
2300 if (me.eq.king .or. .not.out1file)
2301 & open(iout,file=outname,status='unknown')
2303 if (fg_rank.gt.0) then
2304 write (liczba,'(i3.3)') myrank/nfgtasks
2305 write (ll,'(bz,i3.3)') fg_rank
2306 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2311 open(igeom,file=intname,status='unknown',access='append')
2312 open(ipdb,file=pdbname,status='unknown')
2313 open(imol2,file=mol2name,status='unknown')
2314 open(istat,file=statname,status='unknown',access='append')
2316 c1out open(iout,file=outname,status='unknown')
2319 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2320 csa_seed=prefix(:lenpre)//'.CSA.seed'
2321 csa_history=prefix(:lenpre)//'.CSA.history'
2322 csa_bank=prefix(:lenpre)//'.CSA.bank'
2323 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2324 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2325 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2326 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2327 csa_int=prefix(:lenpre)//'.int'
2328 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2329 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2330 csa_in=prefix(:lenpre)//'.CSA.in'
2331 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2334 write (iout,'(80(1h-))')
2335 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2336 write (iout,'(80(1h-))')
2337 write (iout,*) "Input file : ",
2338 & pref_orig(:ilen(pref_orig))//'.inp'
2339 write (iout,*) "Output file : ",
2340 & outname(:ilen(outname))
2342 write (iout,*) "Sidechain potential file : ",
2343 & sidename(:ilen(sidename))
2345 write (iout,*) "SCp potential file : ",
2346 & scpname(:ilen(scpname))
2348 write (iout,*) "Electrostatic potential file : ",
2349 & elename(:ilen(elename))
2350 write (iout,*) "Cumulant coefficient file : ",
2351 & fouriername(:ilen(fouriername))
2352 write (iout,*) "Torsional parameter file : ",
2353 & torname(:ilen(torname))
2354 write (iout,*) "Double torsional parameter file : ",
2355 & tordname(:ilen(tordname))
2356 write (iout,*) "SCCOR parameter file : ",
2357 & sccorname(:ilen(sccorname))
2358 write (iout,*) "Bond & inertia constant file : ",
2359 & bondname(:ilen(bondname))
2360 write (iout,*) "Bending parameter file : ",
2361 & thetname(:ilen(thetname))
2362 write (iout,*) "Rotamer parameter file : ",
2363 & rotname(:ilen(rotname))
2364 write (iout,*) "Thetpdb parameter file : ",
2365 & thetname_pdb(:ilen(thetname_pdb))
2366 write (iout,*) "Threading database : ",
2367 & patname(:ilen(patname))
2369 &write (iout,*)" DIRTMP : ",
2371 write (iout,'(80(1h-))')
2375 c----------------------------------------------------------------------------
2376 subroutine card_concat(card)
2377 implicit real*8 (a-h,o-z)
2378 include 'DIMENSIONS'
2379 include 'COMMON.IOUNITS'
2381 character*80 karta,ucase
2383 read (inp,'(a)') karta
2386 do while (karta(80:80).eq.'&')
2387 card=card(:ilen(card)+1)//karta(:79)
2388 read (inp,'(a)') karta
2391 card=card(:ilen(card)+1)//karta
2394 c----------------------------------------------------------------------------------
2396 implicit real*8 (a-h,o-z)
2397 include 'DIMENSIONS'
2398 include 'COMMON.CHAIN'
2399 include 'COMMON.IOUNITS'
2401 include 'COMMON.CONTROL'
2402 open(irest2,file=rest2name,status='unknown')
2403 read(irest2,*) totT,EK,potE,totE,t_bath
2406 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2409 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2411 if(usampl.or.homol_nset.gt.1) then
2412 read (irest2,*) iset
2417 c---------------------------------------------------------------------------------
2418 subroutine read_fragments
2419 implicit real*8 (a-h,o-z)
2420 include 'DIMENSIONS'
2424 include 'COMMON.SETUP'
2425 include 'COMMON.CHAIN'
2426 include 'COMMON.IOUNITS'
2428 include 'COMMON.CONTROL'
2430 read(inp,*) nset,nfrag,npair,nfrag_back
2431 if(me.eq.king.or..not.out1file)
2432 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2433 & " nfrag_back",nfrag_back
2435 read(inp,*) mset(iset1)
2437 read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
2439 if(me.eq.king.or..not.out1file)
2440 & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
2441 & ifrag(2,i,iset1), qinfrag(i,iset1)
2444 read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
2446 if(me.eq.king.or..not.out1file)
2447 & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
2448 & ipair(2,i,iset1), qinpair(i,iset1)
2451 read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
2452 & wfrag_back(3,i,iset1),
2453 & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
2454 if(me.eq.king.or..not.out1file)
2455 & write(iout,*) "A",i,wfrag_back(1,i,iset1),
2456 & wfrag_back(2,i,iset1),
2457 & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
2458 & ifrag_back(2,i,iset1)
2463 C---------------------------------------------------------------------------
2464 subroutine read_afminp
2465 implicit real*8 (a-h,o-z)
2466 include 'DIMENSIONS'
2470 include 'COMMON.SETUP'
2471 include 'COMMON.CONTROL'
2472 include 'COMMON.CHAIN'
2473 include 'COMMON.IOUNITS'
2474 include 'COMMON.SBRIDGE'
2475 character*320 afmcard
2476 c print *, "wchodze"
2477 call card_concat(afmcard)
2478 call readi(afmcard,"BEG",afmbeg,0)
2479 call readi(afmcard,"END",afmend,0)
2480 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2481 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2482 print *,'FORCE=' ,forceAFMconst
2483 CCCC NOW PROPERTIES FOR AFM
2486 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2488 distafminit=dsqrt(distafminit)
2489 print *,'initdist',distafminit
2493 c-------------------------------------------------------------------------------
2494 subroutine read_saxs_constr
2495 implicit real*8 (a-h,o-z)
2496 include 'DIMENSIONS'
2500 include 'COMMON.SETUP'
2501 include 'COMMON.CONTROL'
2502 include 'COMMON.CHAIN'
2503 include 'COMMON.IOUNITS'
2504 include 'COMMON.SBRIDGE'
2505 double precision cm(3)
2507 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2509 if (saxs_mode.eq.0) then
2510 c SAXS distance distribution
2512 read(inp,*) distsaxs(i),Psaxs(i)
2516 Cnorm = Cnorm + Psaxs(i)
2518 write (iout,*) "Cnorm",Cnorm
2520 Psaxs(i)=Psaxs(i)/Cnorm
2522 write (iout,*) "Normalized distance distribution from SAXS"
2524 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2528 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2530 write (iout,*) "Wsaxs0",Wsaxs0
2534 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2541 cm(j)=cm(j)+Csaxs(j,i)
2549 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2552 write (iout,*) "SAXS sphere coordinates"
2554 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2559 c-------------------------------------------------------------------------------
2560 subroutine read_dist_constr
2561 implicit real*8 (a-h,o-z)
2562 include 'DIMENSIONS'
2566 include 'COMMON.SETUP'
2567 include 'COMMON.CONTROL'
2568 include 'COMMON.CHAIN'
2569 include 'COMMON.IOUNITS'
2570 include 'COMMON.SBRIDGE'
2571 integer ifrag_(2,100),ipair_(2,1000)
2572 double precision wfrag_(100),wpair_(1000)
2573 character*5000 controlcard
2574 logical normalize,next
2576 double precision xlink(4,0:4) /
2578 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
2579 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
2580 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
2581 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
2582 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
2583 c print *, "WCHODZE"
2584 write (iout,*) "Calling read_dist_constr"
2585 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2595 call card_concat(controlcard)
2596 next = index(controlcard,"NEXT").gt.0
2597 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2598 write (iout,*) "restr_type",restr_type
2599 call readi(controlcard,"NFRAG",nfrag_,0)
2600 call readi(controlcard,"NFRAG",nfrag_,0)
2601 call readi(controlcard,"NPAIR",npair_,0)
2602 call readi(controlcard,"NDIST",ndist_,0)
2603 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2604 if (restr_type.eq.10)
2605 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2606 if (restr_type.eq.12)
2607 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2608 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2609 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2610 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2611 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2612 normalize = index(controlcard,"NORMALIZE").gt.0
2613 write (iout,*) "WBOLTZD",wboltzd
2614 write (iout,*) "SCAL_PEAK",scal_peak
2615 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2616 write (iout,*) "IFRAG"
2618 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2620 write (iout,*) "IPAIR"
2622 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2624 if (nfrag_.gt.0) write (iout,*)
2625 & "Distance restraints as generated from reference structure"
2627 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2628 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2629 & ifrag_(2,i)=nstart_sup+nsup-1
2630 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2632 if (wfrag_(i).eq.0.0d0) cycle
2633 do j=ifrag_(1,i),ifrag_(2,i)-1
2634 do k=j+1,ifrag_(2,i)
2635 c write (iout,*) "j",j," k",k
2637 if (restr_type.eq.1) then
2643 forcon(nhpb)=wfrag_(i)
2644 else if (constr_dist.eq.2) then
2645 if (ddjk.le.dist_cut) then
2651 forcon(nhpb)=wfrag_(i)
2653 else if (restr_type.eq.3) then
2659 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2662 if (.not.out1file .or. me.eq.king)
2663 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2664 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2666 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2667 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2673 if (wpair_(i).eq.0.0d0) cycle
2681 do j=ifrag_(1,ii),ifrag_(2,ii)
2682 do k=ifrag_(1,jj),ifrag_(2,jj)
2684 if (restr_type.eq.1) then
2690 forcon(nhpb)=wpair_(i)
2691 else if (constr_dist.eq.2) then
2692 if (ddjk.le.dist_cut) then
2698 forcon(nhpb)=wpair_(i)
2700 else if (restr_type.eq.3) then
2706 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2709 if (.not.out1file .or. me.eq.king)
2710 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2711 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2713 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2714 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2721 write (iout,*) "Distance restraints as read from input"
2723 if (restr_type.eq.12) then
2724 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2725 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2726 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2727 & fordepth_peak(nhpb_peak+1),npeak
2728 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2729 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2730 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2731 c & fordepth_peak(nhpb_peak+1),npeak
2732 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2733 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2734 nhpb_peak=nhpb_peak+1
2735 irestr_type_peak(nhpb_peak)=12
2736 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2739 if (.not.out1file .or. me.eq.king)
2740 & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2741 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2742 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2743 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2744 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2746 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2747 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2748 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2749 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2750 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2752 if (ibecarb_peak(nhpb_peak).gt.0) then
2753 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2754 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2756 else if (restr_type.eq.11) then
2757 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2758 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2759 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2760 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2762 irestr_type(nhpb)=11
2764 if (.not.out1file .or. me.eq.king)
2765 & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2766 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2767 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2769 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2770 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2771 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2773 if (ibecarb(nhpb).gt.0) then
2774 ihpb(nhpb)=ihpb(nhpb)+nres
2775 jhpb(nhpb)=jhpb(nhpb)+nres
2777 else if (restr_type.eq.10) then
2778 c Cross-lonk Markov-like potential
2779 call card_concat(controlcard)
2780 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2781 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2783 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2784 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2785 if (index(controlcard,"ZL").gt.0) then
2787 else if (index(controlcard,"ADH").gt.0) then
2789 else if (index(controlcard,"PDH").gt.0) then
2791 else if (index(controlcard,"DSS").gt.0) then
2796 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2797 & xlink(1,link_type))
2798 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2799 & xlink(2,link_type))
2800 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2801 & xlink(3,link_type))
2802 call reada(controlcard,"SIGMA",forcon(nhpb+1),
2803 & xlink(4,link_type))
2804 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2805 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2806 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2807 if (forcon(nhpb+1).le.0.0d0 .or.
2808 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2810 irestr_type(nhpb)=10
2811 if (ibecarb(nhpb).gt.0) then
2812 ihpb(nhpb)=ihpb(nhpb)+nres
2813 jhpb(nhpb)=jhpb(nhpb)+nres
2816 if (.not.out1file .or. me.eq.king)
2817 & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2818 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2819 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2822 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2823 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2824 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2829 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2830 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2831 if (forcon(nhpb+1).gt.0.0d0) then
2833 if (dhpb1(nhpb).eq.0.0d0) then
2838 if (ibecarb(nhpb).gt.0) then
2839 ihpb(nhpb)=ihpb(nhpb)+nres
2840 jhpb(nhpb)=jhpb(nhpb)+nres
2842 if (dhpb(nhpb).eq.0.0d0)
2843 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2846 if (.not.out1file .or. me.eq.king)
2847 & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2848 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2850 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2851 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2854 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2855 C if (forcon(nhpb+1).gt.0.0d0) then
2857 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2865 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
2866 & fordepthmax=fordepth(i)
2869 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
2874 c-------------------------------------------------------------------------------
2876 subroutine read_constr_homology
2878 include 'DIMENSIONS'
2882 include 'COMMON.SETUP'
2883 include 'COMMON.CONTROL'
2884 include 'COMMON.CHAIN'
2885 include 'COMMON.IOUNITS'
2887 include 'COMMON.GEO'
2888 include 'COMMON.INTERACT'
2889 include 'COMMON.NAMES'
2891 c For new homol impl
2893 include 'COMMON.VAR'
2896 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2898 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2899 c & sigma_odl_temp(maxres,maxres,max_template)
2901 character*24 model_ki_dist, model_ki_angle
2902 character*500 controlcard
2903 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
2904 logical lprn /.true./
2909 c FP - Nov. 2014 Temporary specifications for new vars
2911 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
2913 double precision, dimension (max_template,maxres) :: rescore
2914 double precision, dimension (max_template,maxres) :: rescore2
2915 double precision, dimension (max_template,maxres) :: rescore3
2916 character*24 pdbfile,tpl_k_rescore
2917 c -----------------------------------------------------------------
2918 c Reading multiple PDB ref structures and calculation of retraints
2919 c not using pre-computed ones stored in files model_ki_{dist,angle}
2921 c -----------------------------------------------------------------
2924 c Alternative: reading from input
2925 call card_concat(controlcard)
2926 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2927 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2928 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2929 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2930 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2931 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
2932 dist1cut=(index(controlcard,'DIST1CUT').gt.0)
2933 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2934 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
2935 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
2936 if(.not.read2sigma.and.start_from_model) then
2937 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
2938 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
2939 start_from_model=.false.
2941 if(start_from_model .and. (me.eq.king .or. .not. out1file))
2942 & write(iout,*) 'START_FROM_MODELS is ON'
2943 if(start_from_model .and. rest) then
2944 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2945 write(iout,*) 'START_FROM_MODELS is OFF'
2946 write(iout,*) 'remove restart keyword from input'
2949 if (homol_nset.gt.1)then
2950 call card_concat(controlcard)
2951 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2952 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2953 write(iout,*) "iset homology_weight "
2955 write(iout,*) i,waga_homology(i)
2958 iset=mod(kolor,homol_nset)+1
2961 waga_homology(1)=1.0
2964 cd write (iout,*) "nnt",nnt," nct",nct
2971 c write(iout,*) 'nnt=',nnt,'nct=',nct
2974 do k=1,constr_homology
2987 if (read_homol_frag) then
2988 call read_klapaucjusz
2991 do k=1,constr_homology
2993 read(inp,'(a)') pdbfile
2994 if(me.eq.king .or. .not. out1file)
2995 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
2996 & pdbfile(:ilen(pdbfile))
2997 open(ipdbin,file=pdbfile,status='old',err=33)
2999 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3000 & pdbfile(:ilen(pdbfile))
3003 c print *,'Begin reading pdb data'
3005 c Files containing res sim or local scores (former containing sigmas)
3008 write(kic2,'(bz,i2.2)') k
3010 tpl_k_rescore="template"//kic2//".sco"
3013 if (read2sigma) then
3014 call readpdb_template(k)
3019 c Distance restraints
3022 C Copy the coordinates from reference coordinates (?)
3026 c write (iout,*) "c(",j,i,") =",c(j,i)
3030 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3032 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3033 open (ientin,file=tpl_k_rescore,status='old')
3034 if (nnt.gt.1) rescore(k,1)=0.0d0
3035 do irec=nnt,nct ! loop for reading res sim
3036 if (read2sigma) then
3037 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3038 & rescore3_tmp,idomain_tmp
3040 idomain(k,i_tmp)=idomain_tmp
3041 rescore(k,i_tmp)=rescore_tmp
3042 rescore2(k,i_tmp)=rescore2_tmp
3043 rescore3(k,i_tmp)=rescore3_tmp
3044 if (.not. out1file .or. me.eq.king)
3045 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3046 & i_tmp,rescore2_tmp,rescore_tmp,
3047 & rescore3_tmp,idomain_tmp
3050 read (ientin,*,end=1401) rescore_tmp
3052 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3053 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3054 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3059 if (waga_dist.ne.0.0d0) then
3067 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3068 c write (iout,*) k,i,j,distal,dist2_cut
3070 if (dist1cut .and. k.gt.1) then
3072 if (l_homo(1,ii)) then
3078 sigma_odl(k,ii)=sigma_odl(1,ii)
3080 l_homo(k,ii)=.false.
3083 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3084 & .and. distal.le.dist2_cut ) then
3090 c write (iout,*) "k",k
3091 c write (iout,*) "i",i," j",j," constr_homology",
3096 if (read2sigma) then
3099 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3101 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3102 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3103 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3105 if (odl(k,ii).le.dist_cut) then
3106 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3109 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3110 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3112 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3113 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3117 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3120 l_homo(k,ii)=.false.
3128 c Theta, dihedral and SC retraints
3130 if (waga_angle.gt.0.0d0) then
3131 c open (ientin,file=tpl_k_sigma_dih,status='old')
3132 c do irec=1,maxres-3 ! loop for reading sigma_dih
3133 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3134 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3135 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3136 c & sigma_dih(k,i+nnt-1)
3141 if (idomain(k,i).eq.0) then
3145 dih(k,i)=phiref(i) ! right?
3146 c read (ientin,*) sigma_dih(k,i) ! original variant
3147 c write (iout,*) "dih(",k,i,") =",dih(k,i)
3148 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3149 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3150 c & "rescore(",k,i-2,") =",rescore(k,i-2),
3151 c & "rescore(",k,i-3,") =",rescore(k,i-3)
3153 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3154 & rescore(k,i-2)+rescore(k,i-3))/4.0
3155 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3156 c write (iout,*) "Raw sigmas for dihedral angle restraints"
3157 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3158 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3159 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
3160 c Instead of res sim other local measure of b/b str reliability possible
3161 if (sigma_dih(k,i).ne.0)
3162 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3163 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3168 if (waga_theta.gt.0.0d0) then
3169 c open (ientin,file=tpl_k_sigma_theta,status='old')
3170 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3171 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3172 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3173 c & sigma_theta(k,i+nnt-1)
3178 do i = nnt+2,nct ! right? without parallel.
3179 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3180 c do i=ithet_start,ithet_end ! with FG parallel.
3181 if (idomain(k,i).eq.0) then
3182 sigma_theta(k,i)=0.0
3185 thetatpl(k,i)=thetaref(i)
3186 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3187 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3188 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3189 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3190 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3191 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3192 & rescore(k,i-2))/3.0
3193 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3194 if (sigma_theta(k,i).ne.0)
3195 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3197 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3198 c rescore(k,i-2) ! right expression ?
3199 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3203 if (waga_d.gt.0.0d0) then
3204 c open (ientin,file=tpl_k_sigma_d,status='old')
3205 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3206 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3207 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3208 c & sigma_d(k,i+nnt-1)
3212 do i = nnt,nct ! right? without parallel.
3213 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3214 c do i=loc_start,loc_end ! with FG parallel.
3215 if (itype(i).eq.10) cycle
3216 if (idomain(k,i).eq.0 ) then
3223 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3224 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3225 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3226 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3227 sigma_d(k,i)=rescore3(k,i) ! right expression ?
3228 if (sigma_d(k,i).ne.0)
3229 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3231 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3232 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3233 c read (ientin,*) sigma_d(k,i) ! 1st variant
3238 c remove distance restraints not used in any model from the list
3239 c shift data in all arrays
3241 if (waga_dist.ne.0.0d0) then
3247 if (ii_in_use(ii).eq.0.and.liiflag) then
3251 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3252 & .not.liiflag.and.ii.eq.lim_odl) then
3253 if (ii.eq.lim_odl) then
3254 iishift=ii-iistart+1
3259 do ki=iistart,lim_odl-iishift
3260 ires_homo(ki)=ires_homo(ki+iishift)
3261 jres_homo(ki)=jres_homo(ki+iishift)
3262 ii_in_use(ki)=ii_in_use(ki+iishift)
3263 do k=1,constr_homology
3264 odl(k,ki)=odl(k,ki+iishift)
3265 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3266 l_homo(k,ki)=l_homo(k,ki+iishift)
3270 lim_odl=lim_odl-iishift
3276 endif ! .not. klapaucjusz
3278 if (constr_homology.gt.0) call homology_partition
3279 if (constr_homology.gt.0) call init_int_table
3280 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3281 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3285 if (.not.lprn) return
3286 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3287 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3288 write (iout,*) "Distance restraints from templates"
3290 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3291 & ii,ires_homo(ii),jres_homo(ii),
3292 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3293 & ki=1,constr_homology)
3295 write (iout,*) "Dihedral angle restraints from templates"
3297 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3298 & (rad2deg*dih(ki,i),
3299 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3301 write (iout,*) "Virtual-bond angle restraints from templates"
3303 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3304 & (rad2deg*thetatpl(ki,i),
3305 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3307 write (iout,*) "SC restraints from templates"
3309 write(iout,'(i5,100(4f8.2,4x))') i,
3310 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3311 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3314 c -----------------------------------------------------------------
3317 c----------------------------------------------------------------------
3320 subroutine flush(iu)
3325 subroutine flush(iu)
3330 c------------------------------------------------------------------------------
3331 subroutine copy_to_tmp(source)
3332 include "DIMENSIONS"
3333 include "COMMON.IOUNITS"
3334 character*(*) source
3335 character* 256 tmpfile
3339 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3340 inquire(file=tmpfile,exist=ex)
3342 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3343 & " to temporary directory..."
3344 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3345 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3349 c------------------------------------------------------------------------------
3350 subroutine move_from_tmp(source)
3351 include "DIMENSIONS"
3352 include "COMMON.IOUNITS"
3353 character*(*) source
3356 write (*,*) "Moving ",source(:ilen(source)),
3357 & " from temporary directory to working directory"
3358 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3359 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3362 c------------------------------------------------------------------------------
3363 subroutine random_init(seed)
3365 C Initialize random number generator
3367 implicit real*8 (a-h,o-z)
3368 include 'DIMENSIONS'
3374 logical OKRandom, prng_restart
3376 integer iseed_array(4)
3378 include 'COMMON.IOUNITS'
3379 include 'COMMON.TIME1'
3380 include 'COMMON.THREAD'
3381 include 'COMMON.SBRIDGE'
3382 include 'COMMON.CONTROL'
3383 include 'COMMON.MCM'
3384 include 'COMMON.MAP'
3385 include 'COMMON.HEADER'
3386 include 'COMMON.CSA'
3387 include 'COMMON.CHAIN'
3388 include 'COMMON.MUCA'
3390 include 'COMMON.FFIELD'
3391 include 'COMMON.SETUP'
3392 iseed=-dint(dabs(seed))
3393 if (iseed.eq.0) then
3394 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3395 & 'Random seed undefined. The program will stop.'
3396 write (*,'(/80(1h*)/20x,a/80(1h*))')
3397 & 'Random seed undefined. The program will stop.'
3399 call mpi_finalize(mpi_comm_world,ierr)
3401 stop 'Bad random seed.'
3404 if (fg_rank.eq.0) then
3408 if(me.eq.king .or. .not. out1file)
3409 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3410 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3411 OKRandom = prng_restart(me,iseedi8)
3414 tmp=65536.0d0**(4-i)
3415 iseed_array(i) = dint(seed/tmp)
3416 seed=seed-iseed_array(i)*tmp
3419 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3420 & (iseed_array(i),i=1,4)
3421 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3422 & (iseed_array(i),i=1,4)
3423 OKRandom = prng_restart(me,iseed_array)
3426 c r1 = prng_next(me)
3427 r1=ran_number(0.0D0,1.0D0)
3428 if(me.eq.king .or. .not. out1file)
3429 & write (iout,*) 'ran_num',r1
3430 if (r1.lt.0.0d0) OKRandom=.false.
3432 if (.not.OKRandom) then
3433 write (iout,*) 'PRNG IS NOT WORKING!!!'
3434 print *,'PRNG IS NOT WORKING!!!'
3437 call mpi_abort(mpi_comm_world,error_msg,ierr)
3440 write (iout,*) 'too many processors for parallel prng'
3441 write (*,*) 'too many processors for parallel prng'
3449 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3454 c -----------------------------------------------------------------
3455 subroutine read_klapaucjusz
3457 include 'DIMENSIONS'
3461 include 'COMMON.SETUP'
3462 include 'COMMON.CONTROL'
3463 include 'COMMON.CHAIN'
3464 include 'COMMON.IOUNITS'
3466 include 'COMMON.GEO'
3467 include 'COMMON.INTERACT'
3468 include 'COMMON.NAMES'
3469 character*256 fragfile
3470 integer ninclust(maxclust),inclust(max_template,maxclust),
3471 & nresclust(maxclust),iresclust(maxres,maxclust)
3474 character*24 model_ki_dist, model_ki_angle
3475 character*500 controlcard
3476 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
3477 logical lprn /.true./
3483 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3484 double precision, dimension (max_template,maxres) :: rescore
3485 double precision, dimension (max_template,maxres) :: rescore2
3486 character*24 pdbfile,tpl_k_rescore
3489 c For new homol impl
3491 include 'COMMON.VAR'
3493 call getenv("FRAGFILE",fragfile)
3494 open(ientin,file=fragfile,status="old",err=10)
3495 read(ientin,*) constr_homology,nclust
3501 do k=1,constr_homology
3502 read(ientin,'(a)') pdbfile
3503 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3504 & pdbfile(:ilen(pdbfile))
3505 open(ipdbin,file=pdbfile,status='old',err=33)
3507 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3508 & pdbfile(:ilen(pdbfile))
3512 call readpdb_template(k)
3520 read(ientin,*) ninclust(i),nresclust(i)
3521 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3522 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3525 c Loop over clusters
3528 do ll = 1,ninclust(l)
3536 idomain(k,iresclust(i,l)+1) = 1
3538 idomain(k,iresclust(i,l)) = 1
3542 c Distance restraints
3545 C Copy the coordinates from reference coordinates (?)
3549 c write (iout,*) "c(",j,i,") =",c(j,i)
3552 call int_from_cart(.true.,.false.)
3553 call sc_loc_geom(.false.)
3555 thetaref(i)=theta(i)
3558 if (waga_dist.ne.0.0d0) then
3566 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3567 c write (iout,*) k,i,j,distal,dist2_cut
3569 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3570 & .and. distal.le.dist2_cut ) then
3576 c write (iout,*) "k",k
3577 c write (iout,*) "i",i," j",j," constr_homology",
3582 if (read2sigma) then
3585 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3587 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3588 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3589 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3591 if (odl(k,ii).le.dist_cut) then
3592 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3595 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3596 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3598 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3599 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3603 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3606 c l_homo(k,ii)=.false.
3613 c Theta, dihedral and SC retraints
3615 if (waga_angle.gt.0.0d0) then
3617 if (idomain(k,i).eq.0) then
3618 c sigma_dih(k,i)=0.0
3622 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3623 & rescore(k,i-2)+rescore(k,i-3))/4.0
3624 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3625 c & " sigma_dihed",sigma_dih(k,i)
3626 if (sigma_dih(k,i).ne.0)
3627 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3632 if (waga_theta.gt.0.0d0) then
3634 if (idomain(k,i).eq.0) then
3635 c sigma_theta(k,i)=0.0
3638 thetatpl(k,i)=thetaref(i)
3639 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3640 & rescore(k,i-2))/3.0
3641 if (sigma_theta(k,i).ne.0)
3642 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3646 if (waga_d.gt.0.0d0) then
3648 if (itype(i).eq.10) cycle
3649 if (idomain(k,i).eq.0 ) then
3656 sigma_d(k,i)=rescore(k,i)
3657 if (sigma_d(k,i).ne.0)
3658 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3659 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3665 c remove distance restraints not used in any model from the list
3666 c shift data in all arrays
3668 if (waga_dist.ne.0.0d0) then
3674 if (ii_in_use(ii).eq.0.and.liiflag) then
3678 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3679 & .not.liiflag.and.ii.eq.lim_odl) then
3680 if (ii.eq.lim_odl) then
3681 iishift=ii-iistart+1
3686 do ki=iistart,lim_odl-iishift
3687 ires_homo(ki)=ires_homo(ki+iishift)
3688 jres_homo(ki)=jres_homo(ki+iishift)
3689 ii_in_use(ki)=ii_in_use(ki+iishift)
3690 do k=1,constr_homology
3691 odl(k,ki)=odl(k,ki+iishift)
3692 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3693 l_homo(k,ki)=l_homo(k,ki+iishift)
3697 lim_odl=lim_odl-iishift
3704 10 stop "Error infragment file"
3706 c----------------------------------------------------------------------