8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.SPLITELE'
14 C Read job setup parameters
16 C Read force-field parameters except weights
18 C Read control parameters for energy minimzation if required
19 if (minim) call read_minim
20 C Read MCM control parameters if required
21 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
22 C Read MD control parameters if reqjuired
23 if (modecalc.eq.12) call read_MDpar
24 C Read MREMD control parameters if required
25 if (modecalc.eq.14) then
29 C Read MUCA control parameters if required
30 if (lmuca) call read_muca
31 C Read CSA control parameters if required (from fort.40 if exists
32 C otherwise from general input file)
33 if (modecalc.eq.8) then
34 inquire (file="fort.40",exist=file_exist)
35 if (.not.file_exist) call csaread
37 cfmc if (modecalc.eq.10) call mcmfread
38 c Read energy-term weights and disulfide parameters
40 C Read molecule information, molecule geometry, energy-term weights, and
41 C restraints if requested
43 C Print restraint information
45 if (.not. out1file .or. me.eq.king) then
48 write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
49 & "The following",nhpb-nss,
50 & " distance restraints have been imposed:",
51 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
54 write (iout,'(4i5,2f8.2,3f10.5,2i5)')i-nss,ihpb(i),jhpb(i),
55 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
60 write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
61 & "The following",npeak,
62 & " NMR peak restraints have been imposed:",
63 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
64 & " score"," type"," ipeak"
66 do j=ipeak(1,i),ipeak(2,i)
67 write (iout,'(5i5,2f8.2,2f10.5,i5)')i,j,ihpb_peak(j),
68 & jhpb_peak(j),ibecarb_peak(j),dhpb_peak(j),dhpb1_peak(j),
69 & forcon_peak(j),fordepth_peak(i),irestr_type_peak(j)
72 write (iout,*) "The ipeak array"
74 write (iout,'(3i5)' ) i,ipeak(1,i),ipeak(2,i)
80 c print *,"Processor",myrank," leaves READRTNS"
83 C-------------------------------------------------------------------------------
84 subroutine read_control
92 logical OKRandom, prng_restart
95 include 'COMMON.IOUNITS'
96 include 'COMMON.TIME1'
97 include 'COMMON.THREAD'
98 include 'COMMON.SBRIDGE'
99 include 'COMMON.CONTROL'
100 include 'COMMON.SAXS'
103 include 'COMMON.HEADER'
105 include 'COMMON.CHAIN'
106 include 'COMMON.MUCA'
108 include 'COMMON.FFIELD'
109 include 'COMMON.INTERACT'
110 include 'COMMON.SETUP'
111 include 'COMMON.SPLITELE'
112 include 'COMMON.SHIELD'
115 integer KDIAG,ICORFL,IXDR
116 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
117 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
119 character*320 controlcard
120 double precision seed
125 read (INP,'(a)') titel
126 call card_concat(controlcard)
127 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
128 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
129 call reada(controlcard,'SEED',seed,0.0D0)
130 call random_init(seed)
131 C Set up the time limit (caution! The time must be input in minutes!)
132 read_cart=index(controlcard,'READ_CART').gt.0
133 out_cart=index(controlcard,'OUT_CART').gt.0
134 out_int=index(controlcard,'OUT_INT').gt.0
135 gmatout=index(controlcard,'GMATOUT').gt.0
136 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
137 C this variable with_theta_constr is the variable which allow to read and execute the
138 C constrains on theta angles WITH_THETA_CONSTR is the keyword
139 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
140 write (iout,*) "with_theta_constr ",with_theta_constr
141 call readi(controlcard,'NSAXS',nsaxs,0)
142 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
143 call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
144 call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
145 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
146 & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
147 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
148 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
149 out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
150 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
151 call readi(controlcard,'SYM',symetr,1)
152 call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
153 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
154 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
155 c call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
156 c call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
157 c call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
158 c call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
159 c call reada(controlcard,'DRMS',drms,0.1D0)
160 c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
161 c write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
162 c write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
163 c write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
164 c write (iout,'(a,f10.1)')'DRMS = ',drms
165 cc write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
166 c write (iout,'(a,f10.1)') 'Time limit (min):',timlim
168 call readi(controlcard,'NZ_START',nz_start,0)
169 call readi(controlcard,'NZ_END',nz_end,0)
170 c call readi(controlcard,'IZ_SC',iz_sc,0)
172 safety = 60.0d0*safety
174 call readi(controlcard,"INTER_LIST_UPDATE",imatupdate,100)
175 call reada(controlcard,"T_BATH",t_bath,300.0d0)
176 minim=(index(controlcard,'MINIMIZE').gt.0)
177 dccart=(index(controlcard,'CART').gt.0)
178 overlapsc=(index(controlcard,'OVERLAP').gt.0)
179 overlapsc=.not.overlapsc
180 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
181 searchsc=.not.searchsc
182 sideadd=(index(controlcard,'SIDEADD').gt.0)
183 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
184 mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
185 outpdb=(index(controlcard,'PDBOUT').gt.0)
186 outmol2=(index(controlcard,'MOL2OUT').gt.0)
187 pdbref=(index(controlcard,'PDBREF').gt.0)
188 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
189 indpdb=index(controlcard,'PDBSTART')
190 extconf=(index(controlcard,'EXTCONF').gt.0)
191 AFMlog=(index(controlcard,'AFM'))
192 selfguide=(index(controlcard,'SELFGUIDE'))
193 c print *,'AFMlog',AFMlog,selfguide,"KUPA"
194 call readi(controlcard,'TUBEMOD',tubelog,0)
195 c write (iout,*) TUBElog,"TUBEMODE"
196 call readi(controlcard,'IPRINT',iprint,0)
197 C SHIELD keyword sets if the shielding effect of side-chains is used
198 C 0 denots no shielding is used all peptide are equally despite the
199 C solvent accesible area
200 C 1 the newly introduced function
201 C 2 reseved for further possible developement
202 call readi(controlcard,'SHIELD',shield_mode,0)
203 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
204 write(iout,*) "shield_mode",shield_mode
206 call readi(controlcard,'TORMODE',tor_mode,0)
207 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
208 write(iout,*) "torsional and valence angle mode",tor_mode
209 call readi(controlcard,'MAXGEN',maxgen,10000)
210 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
211 call readi(controlcard,"KDIAG",kdiag,0)
212 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
213 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
214 & write (iout,*) "RESCALE_MODE",rescale_mode
215 split_ene=index(controlcard,'SPLIT_ENE').gt.0
216 if (index(controlcard,'REGULAR').gt.0.0D0) then
217 call reada(controlcard,'WEIDIS',weidis,0.1D0)
222 if (index(controlcard,"CASC").gt.0) then
224 c else if (index(controlcard,"CAONLY").gt.0) then
226 else if (index(controlcard,"SCONLY").gt.0) then
230 c write (iout,*) "ERROR! Must specify CASC CAONLY or SCONLY when",
231 c & " specifying REFSTR, PDBREF or REGULAR."
235 if (index(controlcard,'CHECKGRAD').gt.0) then
237 if (index(controlcard,' CART').gt.0) then
239 elseif (index(controlcard,'CARINT').gt.0) then
244 call reada(controlcard,'DELTA',aincr,1.0d-7)
245 c write (iout,*) "icheckgrad",icheckgrad
246 elseif (index(controlcard,'THREAD').gt.0) then
248 call readi(controlcard,'THREAD',nthread,0)
249 if (nthread.gt.0) then
250 call reada(controlcard,'WEIDIS',weidis,0.1D0)
253 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
254 stop 'Error termination in Read_Control.'
256 else if (index(controlcard,'MCMA').gt.0) then
258 else if (index(controlcard,'MCEE').gt.0) then
260 else if (index(controlcard,'MULTCONF').gt.0) then
262 else if (index(controlcard,'MAP').gt.0) then
264 call readi(controlcard,'MAP',nmap,0)
265 else if (index(controlcard,'CSA').gt.0) then
267 crc else if (index(controlcard,'ZSCORE').gt.0) then
269 crc ZSCORE is rm from UNRES, modecalc=9 is available
272 cfcm else if (index(controlcard,'MCMF').gt.0) then
274 else if (index(controlcard,'SOFTREG').gt.0) then
276 else if (index(controlcard,'CHECK_BOND').gt.0) then
278 else if (index(controlcard,'TEST').gt.0) then
280 else if (index(controlcard,'MD').gt.0) then
282 else if (index(controlcard,'RE ').gt.0) then
286 lmuca=index(controlcard,'MUCA').gt.0
287 call readi(controlcard,'MUCADYN',mucadyn,0)
288 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
289 if (lmuca .and. (me.eq.king .or. .not.out1file ))
291 write (iout,*) 'MUCADYN=',mucadyn
292 write (iout,*) 'MUCASMOOTH=',muca_smooth
295 iscode=index(controlcard,'ONE_LETTER')
296 indphi=index(controlcard,'PHI')
297 indback=index(controlcard,'BACK')
298 iranconf=index(controlcard,'RAND_CONF')
299 i2ndstr=index(controlcard,'USE_SEC_PRED')
300 gradout=index(controlcard,'GRADOUT').gt.0
301 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
302 C DISTCHAINMAX become obsolete for periodic boundry condition
303 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
304 C Reading the dimensions of box in x,y,z coordinates
305 call reada(controlcard,'BOXX',boxxsize,100.0d0)
306 call reada(controlcard,'BOXY',boxysize,100.0d0)
307 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
308 write(iout,*) "Periodic box dimensions",boxxsize,boxysize,boxzsize
309 c Cutoff range for interactions
310 call reada(controlcard,"R_CUT_INT",r_cut_int,25.0d0)
311 call reada(controlcard,"R_CUT_RESPA",r_cut_respa,2.0d0)
312 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
313 write (iout,*) "Cutoff on interactions",r_cut_int
315 & "Cutoff in switching short and long range interactions in RESPA",
317 write (iout,*) "lambda in switch function",rlamb
318 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
319 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
320 if (lipthick.gt.0.0d0) then
321 bordliptop=(boxzsize+lipthick)/2.0
322 bordlipbot=bordliptop-lipthick
324 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
325 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
326 buflipbot=bordlipbot+lipbufthick
327 bufliptop=bordliptop-lipbufthick
328 if ((lipbufthick*2.0d0).gt.lipthick)
329 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
331 write(iout,*) "bordliptop=",bordliptop
332 write(iout,*) "bordlipbot=",bordlipbot
333 write(iout,*) "bufliptop=",bufliptop
334 write(iout,*) "buflipbot=",buflipbot
335 write (iout,*) "SHIELD MODE",shield_mode
336 if (TUBElog.gt.0) then
337 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
338 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
339 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
340 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
341 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
342 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
343 buftubebot=bordtubebot+tubebufthick
344 buftubetop=bordtubetop-tubebufthick
346 if (me.eq.king .or. .not.out1file )
347 & write (iout,*) "DISTCHAINMAX",distchainmax
349 if(me.eq.king.or..not.out1file)
350 & write (iout,'(2a)') diagmeth(kdiag),
351 & ' routine used to diagonalize matrices.'
354 c--------------------------------------------------------------------------
355 subroutine read_REMDpar
361 include 'COMMON.IOUNITS'
362 include 'COMMON.TIME1'
365 include 'COMMON.LANGEVIN'
368 include 'COMMON.LANGEVIN.lang0.5diag'
370 include 'COMMON.LANGEVIN.lang0'
373 include 'COMMON.INTERACT'
374 include 'COMMON.NAMES'
376 include 'COMMON.REMD'
377 include 'COMMON.CONTROL'
378 include 'COMMON.SETUP'
380 character*320 controlcard
381 character*3200 controlcard1
382 integer iremd_m_total,i
384 if(me.eq.king.or..not.out1file)
385 & write (iout,*) "REMD setup"
387 call card_concat(controlcard)
388 call readi(controlcard,"NREP",nrep,3)
389 call readi(controlcard,"NSTEX",nstex,1000)
390 call reada(controlcard,"RETMIN",retmin,10.0d0)
391 call reada(controlcard,"RETMAX",retmax,1000.0d0)
392 mremdsync=(index(controlcard,'SYNC').gt.0)
393 call readi(controlcard,"NSYN",i_sync_step,100)
394 restart1file=(index(controlcard,'REST1FILE').gt.0)
395 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
396 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
397 if(max_cache_traj_use.gt.max_cache_traj)
398 & max_cache_traj_use=max_cache_traj
399 if(me.eq.king.or..not.out1file) then
400 cd if (traj1file) then
401 crc caching is in testing - NTWX is not ignored
402 cd write (iout,*) "NTWX value is ignored"
403 cd write (iout,*) " trajectory is stored to one file by master"
404 cd write (iout,*) " before exchange at NSTEX intervals"
406 write (iout,*) "NREP= ",nrep
407 write (iout,*) "NSTEX= ",nstex
408 write (iout,*) "SYNC= ",mremdsync
409 write (iout,*) "NSYN= ",i_sync_step
410 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
413 if (index(controlcard,'TLIST').gt.0) then
415 call card_concat(controlcard1)
416 read(controlcard1,*) (remd_t(i),i=1,nrep)
417 if(me.eq.king.or..not.out1file)
418 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
421 if (index(controlcard,'MLIST').gt.0) then
423 call card_concat(controlcard1)
424 read(controlcard1,*) (remd_m(i),i=1,nrep)
425 if(me.eq.king.or..not.out1file) then
426 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
429 iremd_m_total=iremd_m_total+remd_m(i)
431 write (iout,*) 'Total number of replicas ',iremd_m_total
436 & "Adaptive (PMF-biased) umbrella sampling will be run"
439 if(me.eq.king.or..not.out1file)
440 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
443 c--------------------------------------------------------------------------
444 subroutine read_MDpar
450 include 'COMMON.IOUNITS'
451 include 'COMMON.TIME1'
453 include 'COMMON.QRESTR'
455 include 'COMMON.LANGEVIN'
458 include 'COMMON.LANGEVIN.lang0.5diag'
460 include 'COMMON.LANGEVIN.lang0'
463 include 'COMMON.INTERACT'
464 include 'COMMON.NAMES'
466 include 'COMMON.SETUP'
467 include 'COMMON.CONTROL'
468 include 'COMMON.SPLITELE'
469 include 'COMMON.FFIELD'
471 character*320 controlcard
475 call card_concat(controlcard)
476 call readi(controlcard,"NSTEP",n_timestep,1000000)
477 call readi(controlcard,"NTWE",ntwe,100)
478 call readi(controlcard,"NTWX",ntwx,1000)
479 call reada(controlcard,"DT",d_time,1.0d-1)
480 call reada(controlcard,"DVMAX",dvmax,2.0d1)
481 call reada(controlcard,"DAMAX",damax,1.0d1)
482 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
483 call readi(controlcard,"LANG",lang,0)
484 RESPA = index(controlcard,"RESPA") .gt. 0
485 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
486 ntime_split0=ntime_split
487 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
488 ntime_split0=ntime_split
489 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
490 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
491 rest = index(controlcard,"REST").gt.0
492 tbf = index(controlcard,"TBF").gt.0
493 usampl = index(controlcard,"USAMPL").gt.0
494 scale_umb = index(controlcard,"SCALE_UMB").gt.0
495 adaptive = index(controlcard,"ADAPTIVE").gt.0
496 mdpdb = index(controlcard,"MDPDB").gt.0
497 call reada(controlcard,"T_BATH",t_bath,300.0d0)
498 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
499 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
500 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
501 if (count_reset_moment.eq.0) count_reset_moment=1000000000
502 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
503 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
504 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
505 if (count_reset_vel.eq.0) count_reset_vel=1000000000
506 large = index(controlcard,"LARGE").gt.0
507 print_compon = index(controlcard,"PRINT_COMPON").gt.0
508 rattle = index(controlcard,"RATTLE").gt.0
509 preminim = index(controlcard,"PREMINIM").gt.0
510 if (iranconf.gt.0 .or. indpdb.gt.0 .or. start_from_model) then
514 write (iout,*) "PREMINIM ",preminim
516 if (index(controlcard,'CART').gt.0) dccart=.true.
517 write (iout,*) "dccart ",dccart
518 write (iout,*) "read_minim ",index(controlcard,'READ_MINIM_PAR')
519 if (index(controlcard,'READ_MINIM_PAR').gt.0) call read_minim
521 c if performing umbrella sampling, fragments constrained are read from the fragment file
524 write (iout,*) "Umbrella sampling will be run"
525 if (scale_umb.and.adaptive) then
526 write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
527 write (iout,*) "Select one of those and re-run the job."
530 if (scale_umb) write (iout,*)
531 &"Umbrella-restraint force constants will be scaled by temperature"
535 if(me.eq.king.or..not.out1file) then
537 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
539 write (iout,'(a)') "The units are:"
540 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
541 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
542 & " acceleration: angstrom/(48.9 fs)**2"
543 write (iout,'(a)') "energy: kcal/mol, temperature: K"
545 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
546 write (iout,'(a60,f10.5,a)')
547 & "Initial time step of numerical integration:",d_time,
549 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
550 write (iout,'(a60,f10.5,a)') "Cutoff on interactions",r_cut_int,
552 write(iout,'(a60,i5)')"Frequency of updating interaction list",
555 write (iout,'(2a,i4,a)')
556 & "A-MTS algorithm used; initial time step for fast-varying",
557 & " short-range forces split into",ntime_split," steps."
558 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
559 & r_cut_respa," lambda",rlamb
561 write (iout,'(2a,f10.5)')
562 & "Maximum acceleration threshold to reduce the time step",
563 & "/increase split number:",damax
564 write (iout,'(2a,f10.5)')
565 & "Maximum predicted energy drift to reduce the timestep",
566 & "/increase split number:",edriftmax
567 write (iout,'(a60,f10.5)')
568 & "Maximum velocity threshold to reduce velocities:",dvmax
569 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
570 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
571 if (rattle) write (iout,'(a60)')
572 & "Rattle algorithm used to constrain the virtual bonds"
576 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
577 call reada(controlcard,"RWAT",rwat,1.4d0)
578 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
579 surfarea=index(controlcard,"SURFAREA").gt.0
580 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
581 if(me.eq.king.or..not.out1file)then
582 write (iout,'(/a,$)') "Langevin dynamics calculation"
585 & " with direct integration of Langevin equations"
586 else if (lang.eq.2) then
587 write (iout,'(a/)') " with TINKER stochasic MD integrator"
588 else if (lang.eq.3) then
589 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
590 else if (lang.eq.4) then
591 write (iout,'(a/)') " in overdamped mode"
593 write (iout,'(//a,i5)')
594 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
597 write (iout,'(a60,f10.5)') "Temperature:",t_bath
598 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
599 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
600 write (iout,'(a60,f10.5)')
601 & "Scaling factor of the friction forces:",scal_fric
602 if (surfarea) write (iout,'(2a,i10,a)')
603 & "Friction coefficients will be scaled by solvent-accessible",
604 & " surface area every",reset_fricmat," steps."
606 c Calculate friction coefficients and bounds of stochastic forces
607 eta=6*pi*cPoise*etawat
608 if(me.eq.king.or..not.out1file)
609 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
611 gamp=scal_fric*(pstok+rwat)*eta
612 stdfp=dsqrt(2*Rb*t_bath/d_time)
614 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
615 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
617 if(me.eq.king.or..not.out1file)then
618 write (iout,'(/2a/)')
619 & "Radii of site types and friction coefficients and std's of",
620 & " stochastic forces of fully exposed sites"
621 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
623 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
624 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
628 if(me.eq.king.or..not.out1file)then
629 write (iout,'(a)') "Berendsen bath calculation"
630 write (iout,'(a60,f10.5)') "Temperature:",t_bath
631 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
633 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
634 & count_reset_moment," steps"
636 & write (iout,'(a,i10,a)')
637 & "Velocities will be reset at random every",count_reset_vel,
641 if(me.eq.king.or..not.out1file)
642 & write (iout,'(a31)') "Microcanonical mode calculation"
644 if(me.eq.king.or..not.out1file)then
645 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
647 write(iout,*) "MD running with constraints."
648 write(iout,*) "Equilibration time ", eq_time, " mtus."
649 write(iout,*) "Constraining ", nfrag," fragments."
650 write(iout,*) "Length of each fragment, weight and q0:"
652 write (iout,*) "Set of restraints #",iset
654 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
655 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
657 write(iout,*) "constraints between ", npair, "fragments."
658 write(iout,*) "constraint pairs, weights and q0:"
660 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
661 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
663 write(iout,*) "angle constraints within ", nfrag_back,
664 & "backbone fragments."
666 write(iout,*) "fragment, weights, q0:"
668 write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
669 & ifrag_back(2,i,iset),
670 & wfrag_back(1,i,iset),qin_back(1,i,iset),
671 & wfrag_back(2,i,iset),qin_back(2,i,iset),
672 & wfrag_back(3,i,iset),qin_back(3,i,iset)
675 write(iout,*) "fragment, weights:"
677 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
678 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
679 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
683 iset=mod(kolor,nset)+1
686 if(me.eq.king.or..not.out1file)
687 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
690 c------------------------------------------------------------------------------
693 C Read molecular data.
699 integer error_msg,ierror,ierr,ierrcode
701 include 'COMMON.IOUNITS'
704 include 'COMMON.INTERACT'
705 include 'COMMON.LOCAL'
706 include 'COMMON.NAMES'
707 include 'COMMON.CHAIN'
708 include 'COMMON.FFIELD'
709 include 'COMMON.SBRIDGE'
710 include 'COMMON.HEADER'
711 include 'COMMON.CONTROL'
712 include 'COMMON.SAXS'
713 include 'COMMON.DBASE'
714 include 'COMMON.THREAD'
715 include 'COMMON.CONTACTS'
716 include 'COMMON.TORCNSTR'
717 include 'COMMON.TIME1'
718 include 'COMMON.BOUNDS'
720 include 'COMMON.SETUP'
721 include 'COMMON.SHIELD'
722 character*4 sequence(maxres)
724 double precision x(maxvar)
725 character*256 pdbfile
726 character*400 weightcard
727 character*80 weightcard_t,ucase
728 integer itype_pdb(maxres)
729 common /pizda/ itype_pdb
730 logical seq_comp,fail
731 double precision energia(0:n_ene)
732 double precision secprob(3,maxdih_constr)
734 double precision phihel,phibet,sigmahel,sigmabet
735 integer iti,nsi,maxsi
739 integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2
740 double precision sumv
742 C Read PDB structure if applicable
744 if (indpdb.gt.0 .or. pdbref) then
745 read(inp,'(a)') pdbfile
746 if(me.eq.king.or..not.out1file)
747 & write (iout,'(2a)') 'PDB data will be read from file ',
748 & pdbfile(:ilen(pdbfile))
749 open(ipdbin,file=pdbfile,status='old',err=33)
751 33 write (iout,'(a)') 'Error opening PDB file.'
762 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
763 & (crefjlee(j,i+nres),j=1,3)
766 if(me.eq.king.or..not.out1file)
767 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
768 & ' nstart_sup=',nstart_sup
770 itype_pdb(i)=itype(i)
774 nct=nstart_sup+nsup-1
775 call contact(.false.,ncont_ref,icont_ref,co)
778 if(me.eq.king.or..not.out1file)
779 & write(iout,*)'Adding sidechains'
783 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
786 do while (fail.and.nsi.le.maxsi)
787 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
790 if(fail) write(iout,*)'Adding sidechain failed for res ',
791 & i,' after ',nsi,' trials'
796 if (indpdb.eq.0) then
797 C Read sequence if not taken from the pdb file.
799 c print *,'nres=',nres
800 if (iscode.gt.0) then
801 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
803 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
805 C Convert sequence to numeric code
807 itype(i)=rescode(i,sequence(i),iscode)
811 c print '(20i4)',(itype(i),i=1,nres)
814 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
816 if (itype(i).eq.ntyp1) then
820 else if (iabs(itype(i+1)).ne.20) then
822 else if (iabs(itype(i)).ne.20) then
829 if(me.eq.king.or..not.out1file)then
830 write (iout,*) "ITEL"
832 write (iout,*) i,itype(i),itel(i)
834 print *,'Call Read_Bridge.'
838 cd print *,'NNT=',NNT,' NCT=',NCT
839 call seq2chains(nres,itype,nchain,chain_length,chain_border,
842 chain_border1(2,1)=chain_border(2,1)+1
844 chain_border1(1,i)=chain_border(1,i)-1
845 chain_border1(2,i)=chain_border(2,i)+1
847 chain_border1(1,nchain)=chain_border(1,nchain)-1
848 chain_border1(2,nchain)=nres
849 write(iout,*) "nres",nres," nchain",nchain
851 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
852 & chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
854 call chain_symmetry(nchain,nres,itype,chain_border,
855 & chain_length,npermchain,tabpermchain)
857 c write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
860 write(iout,*) "residue permutations"
862 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
865 if (itype(1).eq.ntyp1) nnt=2
866 if (itype(nres).eq.ntyp1) nct=nct-1
867 write (iout,*) "nnt",nnt," nct",nct
870 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
871 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
873 print*, 'init_dfa_vars finished!'
875 print*, 'read_dfa_info finished!'
879 if(me.eq.king.or..not.out1file)
880 & write (iout,'(a,i3)') 'nsup=',nsup
882 if (nsup.le.(nct-nnt+1)) then
883 do i=0,nct-nnt+1-nsup
884 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
890 & 'Error - sequences to be superposed do not match.'
893 do i=0,nsup-(nct-nnt+1)
894 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
896 nstart_sup=nstart_sup+i
902 & 'Error - sequences to be superposed do not match.'
905 if (nsup.eq.0) nsup=nct-nnt
906 if (nstart_sup.eq.0) nstart_sup=nnt
907 if (nstart_seq.eq.0) nstart_seq=nnt
908 if(me.eq.king.or..not.out1file)
909 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
910 & ' nstart_seq=',nstart_seq
913 C 8/13/98 Set limits to generating the dihedral angles
918 read (inp,*) ndih_constr
919 write (iout,*) "ndih_constr",ndih_constr
920 if (ndih_constr.gt.0) then
923 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
926 if(me.eq.king.or..not.out1file)then
929 & 'There are',ndih_constr,' restraints on gamma angles.'
931 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
938 phi0(i)=deg2rad*phi0(i)
939 drange(i)=deg2rad*drange(i)
941 C if(me.eq.king.or..not.out1file)
942 C & write (iout,*) 'FTORS',ftors
945 phibound(1,ii) = phi0(i)-drange(i)
946 phibound(2,ii) = phi0(i)+drange(i)
949 if (me.eq.king .or. .not.out1file) then
951 write (iout,'(a)') 'Boundaries in gamma angle sampling:'
953 write (iout,'(a3,i5,2f10.1)')
954 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
959 else if (ndih_constr.lt.0) then
961 call card_concat(weightcard)
962 call reada(weightcard,"PHIHEL",phihel,50.0D0)
963 call reada(weightcard,"PHIBET",phibet,180.0D0)
964 call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
965 call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
966 call reada(weightcard,"WDIHC",wdihc,0.591D0)
967 write (iout,*) "Weight of dihedral angle restraints",wdihc
968 read(inp,'(9x,3f7.3)')
969 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
970 write (iout,*) "The secprob array"
972 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
976 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
977 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
978 ndih_constr=ndih_constr+1
979 idih_constr(ndih_constr)=i
982 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
983 sumv=sumv+vpsipred(j,ndih_constr)
986 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
988 phibound(1,ndih_constr)=phihel*deg2rad
989 phibound(2,ndih_constr)=phibet*deg2rad
990 sdihed(1,ndih_constr)=sigmahel*deg2rad
991 sdihed(2,ndih_constr)=sigmabet*deg2rad
995 if(me.eq.king.or..not.out1file)then
998 & 'There are',ndih_constr,
999 & ' bimodal restraints on gamma angles.'
1001 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
1002 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
1003 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
1004 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
1005 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
1006 & (vpsipred(j,i),j=1,3)
1013 C first setting the theta boundaries to 0 to pi
1014 C this mean that there is no energy penalty for any angle occuring this can be applied
1015 C for generate random conformation but is not implemented in this way
1018 C thetabound(2,i)=pi
1020 C begin reading theta constrains this is quartic constrains allowing to
1021 C have smooth second derivative
1022 if (with_theta_constr) then
1023 C with_theta_constr is keyword allowing for occurance of theta constrains
1024 read (inp,*) ntheta_constr
1025 C ntheta_constr is the number of theta constrains
1026 if (ntheta_constr.gt.0) then
1027 C read (inp,*) ftors
1028 read (inp,*) (itheta_constr(i),theta_constr0(i),
1029 & theta_drange(i),for_thet_constr(i),
1030 & i=1,ntheta_constr)
1031 C the above code reads from 1 to ntheta_constr
1032 C itheta_constr(i) residue i for which is theta_constr
1033 C theta_constr0 the global minimum value
1034 C theta_drange is range for which there is no energy penalty
1035 C for_thet_constr is the force constant for quartic energy penalty
1037 if(me.eq.king.or..not.out1file)then
1039 & 'There are',ntheta_constr,' constraints on phi angles.'
1040 do i=1,ntheta_constr
1041 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1043 & for_thet_constr(i)
1046 do i=1,ntheta_constr
1047 theta_constr0(i)=deg2rad*theta_constr0(i)
1048 theta_drange(i)=deg2rad*theta_drange(i)
1050 C if(me.eq.king.or..not.out1file)
1051 C & write (iout,*) 'FTORS',ftors
1052 C do i=1,ntheta_constr
1053 C ii = itheta_constr(i)
1054 C thetabound(1,ii) = phi0(i)-drange(i)
1055 C thetabound(2,ii) = phi0(i)+drange(i)
1057 endif ! ntheta_constr.gt.0
1058 endif! with_theta_constr
1060 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1061 C write (iout,*) "with_dihed_constr ",with_dihed_constr
1062 c--- Zscore rms -------
1063 if (nz_start.eq.0) nz_start=nnt
1064 if (nz_end.eq.0 .and. nsup.gt.0) then
1066 else if (nz_end.eq.0) then
1069 if(me.eq.king.or..not.out1file)then
1070 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1071 write (iout,*) 'IZ_SC=',iz_sc
1073 c----------------------
1076 if (.not.pdbref) then
1077 call read_angles(inp,*38)
1080 38 write (iout,'(a)') 'Error reading reference structure.'
1082 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1083 stop 'Error reading reference structure'
1085 39 call chainbuild_extconf
1087 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1096 call contact(.true.,ncont_ref,icont_ref,co)
1100 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1102 if (constr_dist.gt.0) call read_dist_constr
1103 write (iout,*) "After read_dist_constr nhpb",nhpb
1104 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1106 call NMRpeak_partition
1107 if(me.eq.king.or..not.out1file)write (iout,*) 'Contact order:',co
1109 if(me.eq.king.or..not.out1file)
1110 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1113 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1115 if(me.eq.king.or..not.out1file)
1116 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1117 & icont_ref(1,i),' ',
1118 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1121 write (iout,*) "calling read_saxs_consrtr",nsaxs
1122 if (nsaxs.gt.0) call read_saxs_constr
1124 if (constr_homology.gt.0) then
1125 call read_constr_homology
1126 if (indpdb.gt.0 .or. pdbref) then
1129 c(j,i)=crefjlee(j,i)
1130 cref(j,i)=crefjlee(j,i)
1135 write (iout,*) "sc_loc_geom: Array C"
1137 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1138 & (c(j,i+nres),j=1,3)
1140 write (iout,*) "Array Cref"
1142 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1143 & (cref(j,i+nres),j=1,3)
1146 call int_from_cart1(.false.)
1147 call sc_loc_geom(.false.)
1149 thetaref(i)=theta(i)
1154 dc(j,i)=c(j,i+1)-c(j,i)
1155 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1160 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1161 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1170 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1171 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1172 & modecalc.ne.10) then
1173 C If input structure hasn't been supplied from the PDB file read or generate
1175 if (iranconf.eq.0 .and. .not. extconf) then
1176 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1177 & write (iout,'(a)') 'Initial geometry will be read in.'
1179 read(inp,'(8f10.5)',end=36,err=36)
1180 & ((c(l,k),l=1,3),k=1,nres),
1181 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1182 write (iout,*) "Exit READ_CART"
1183 c write (iout,'(8f10.5)')
1184 c & ((c(l,k),l=1,3),k=1,nres),
1185 c & ((c(l,k+nres),l=1,3),k=nnt,nct)
1189 dc(j,i)=c(j,i+1)-c(j,i)
1190 c dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1194 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1196 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1197 c dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1202 c dc_norm(j,i+nres)=0.0d0
1206 call int_from_cart1(.true.)
1207 write (iout,*) "Finish INT_TO_CART"
1208 c write (iout,*) "DC"
1210 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1211 c & (dc(j,i+nres),j=1,3)
1217 call read_angles(inp,*36)
1219 call chainbuild_extconf
1222 36 write (iout,'(a)') 'Error reading angle file.'
1224 call mpi_finalize( MPI_COMM_WORLD,IERR )
1226 stop 'Error reading angle file.'
1228 else if (extconf) then
1229 if (me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1230 & write (iout,'(a)') 'Extended chain initial geometry.'
1232 theta(i)=90d0*deg2rad
1235 phi(i)=180d0*deg2rad
1238 alph(i)=110d0*deg2rad
1241 omeg(i)=-120d0*deg2rad
1242 if (itype(i).le.0) omeg(i)=-omeg(i)
1245 call chainbuild_extconf
1247 if(me.eq.king.or..not.out1file)
1248 & write (iout,'(a)') 'Random-generated initial geometry.'
1251 if (me.eq.king .or. fg_rank.eq.0 .and. (
1252 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1256 call gen_rand_conf(itmp,*30)
1258 30 write (iout,*) 'Failed to generate random conformation',
1259 & ', itrial=',itrial
1260 write (*,*) 'Processor:',me,
1261 & ' Failed to generate random conformation',
1271 write (iout,'(a,i3,a)') 'Processor:',me,
1272 & ' error in generating random conformation.'
1273 write (*,'(a,i3,a)') 'Processor:',me,
1274 & ' error in generating random conformation.'
1277 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1282 & ' error in generating random conformation.'
1287 elseif (modecalc.eq.4) then
1288 read (inp,'(a)') intinname
1289 open (intin,file=intinname,status='old',err=333)
1290 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1291 & write (iout,'(a)') 'intinname',intinname
1292 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1294 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1296 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1298 stop 'Error opening angle file.'
1302 C Generate distance constraints, if the PDB structure is to be regularized.
1303 if (nthread.gt.0) then
1304 call read_threadbase
1307 if (me.eq.king .or. .not. out1file)
1309 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1310 write (iout,'(/a,i3,a)')
1311 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1312 write (iout,'(20i4)') (iss(i),i=1,ns)
1314 write(iout,*)"Running with dynamic disulfide-bond formation"
1316 write (iout,'(/a/)') 'Pre-formed links are:'
1322 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1323 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1329 if (ns.gt.0.and.dyn_ss) then
1333 forcon(i-nss)=forcon(i)
1340 dyn_ss_mask(iss(i))=.true.
1345 c write (iout,*) "DC"
1347 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1348 c & (dc(j,i+nres),j=1,3)
1350 if (i2ndstr.gt.0) call secstrp2dihc
1351 c call geom_to_var(nvar,x)
1352 c call etotal(energia(0))
1353 c call enerprint(energia(0))
1354 c call briefout(0,etot)
1356 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1357 cd write (iout,'(a)') 'Variable list:'
1358 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1360 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1361 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1362 & 'Processor',myrank,': end reading molecular data.'
1367 c--------------------------------------------------------------------------
1368 logical function seq_comp(itypea,itypeb,length)
1370 integer length,itypea(length),itypeb(length)
1373 if (itypea(i).ne.itypeb(i)) then
1381 c-----------------------------------------------------------------------------
1382 subroutine read_bridge
1383 C Read information about disulfide bridges.
1385 include 'DIMENSIONS'
1390 include 'COMMON.IOUNITS'
1391 include 'COMMON.GEO'
1392 include 'COMMON.VAR'
1393 include 'COMMON.INTERACT'
1394 include 'COMMON.LOCAL'
1395 include 'COMMON.NAMES'
1396 include 'COMMON.CHAIN'
1397 include 'COMMON.FFIELD'
1398 include 'COMMON.SBRIDGE'
1399 include 'COMMON.HEADER'
1400 include 'COMMON.CONTROL'
1401 include 'COMMON.DBASE'
1402 include 'COMMON.THREAD'
1403 include 'COMMON.TIME1'
1404 include 'COMMON.SETUP'
1406 C Read bridging residues.
1407 read (inp,*) ns,(iss(i),i=1,ns)
1409 if(me.eq.king.or..not.out1file)
1410 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1411 C Check whether the specified bridging residues are cystines.
1413 if (itype(iss(i)).ne.1) then
1414 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1415 & 'Do you REALLY think that the residue ',
1416 & restyp(itype(iss(i))),i,
1417 & ' can form a disulfide bridge?!!!'
1418 write (*,'(2a,i3,a)')
1419 & 'Do you REALLY think that the residue ',
1420 & restyp(itype(iss(i))),i,
1421 & ' can form a disulfide bridge?!!!'
1423 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1428 C Read preformed bridges.
1430 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1432 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1435 C Check if the residues involved in bridges are in the specified list of
1436 C bridging residues.
1439 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1440 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1441 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1442 & ' contains residues present in other pairs.'
1443 write (*,'(a,i3,a)') 'Disulfide pair',i,
1444 & ' contains residues present in other pairs.'
1446 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1452 if (ihpb(i).eq.iss(j)) goto 10
1454 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1457 if (jhpb(i).eq.iss(j)) goto 20
1459 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1465 ihpb(i)=ihpb(i)+nres
1466 jhpb(i)=jhpb(i)+nres
1472 c----------------------------------------------------------------------------
1473 subroutine read_x(kanal,*)
1475 include 'DIMENSIONS'
1476 include 'COMMON.GEO'
1477 include 'COMMON.VAR'
1478 include 'COMMON.CHAIN'
1479 include 'COMMON.IOUNITS'
1480 include 'COMMON.CONTROL'
1481 include 'COMMON.LOCAL'
1482 include 'COMMON.INTERACT'
1483 integer i,j,k,l,kanal
1484 c Read coordinates from input
1486 read(kanal,'(8f10.5)',end=10,err=10)
1487 & ((c(l,k),l=1,3),k=1,nres),
1488 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1491 c(j,2*nres)=c(j,nres)
1493 call int_from_cart1(.false.)
1496 dc(j,i)=c(j,i+1)-c(j,i)
1497 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1501 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1503 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1504 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1512 c----------------------------------------------------------------------------
1513 subroutine read_threadbase
1515 include 'DIMENSIONS'
1516 include 'COMMON.IOUNITS'
1517 include 'COMMON.GEO'
1518 include 'COMMON.VAR'
1519 include 'COMMON.INTERACT'
1520 include 'COMMON.LOCAL'
1521 include 'COMMON.NAMES'
1522 include 'COMMON.CHAIN'
1523 include 'COMMON.FFIELD'
1524 include 'COMMON.SBRIDGE'
1525 include 'COMMON.HEADER'
1526 include 'COMMON.CONTROL'
1527 include 'COMMON.DBASE'
1528 include 'COMMON.THREAD'
1529 include 'COMMON.TIME1'
1531 double precision dist
1532 C Read pattern database for threading.
1533 read (icbase,*) nseq
1535 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1536 & nres_base(2,i),nres_base(3,i)
1537 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1539 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1540 c & nres_base(2,i),nres_base(3,i)
1541 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1545 if (weidis.eq.0.0D0) weidis=0.1D0
1554 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1555 write (iout,'(a,i5)') 'nexcl: ',nexcl
1556 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1559 c------------------------------------------------------------------------------
1560 subroutine setup_var
1562 include 'DIMENSIONS'
1563 include 'COMMON.IOUNITS'
1564 include 'COMMON.GEO'
1565 include 'COMMON.VAR'
1566 include 'COMMON.INTERACT'
1567 include 'COMMON.LOCAL'
1568 include 'COMMON.NAMES'
1569 include 'COMMON.CHAIN'
1570 include 'COMMON.FFIELD'
1571 include 'COMMON.SBRIDGE'
1572 include 'COMMON.HEADER'
1573 include 'COMMON.CONTROL'
1574 include 'COMMON.DBASE'
1575 include 'COMMON.THREAD'
1576 include 'COMMON.TIME1'
1578 C Set up variable list.
1584 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1586 ialph(i,1)=nvar+nside
1590 if (indphi.gt.0) then
1592 else if (indback.gt.0) then
1599 c----------------------------------------------------------------------------
1600 subroutine gen_dist_constr
1601 C Generate CA distance constraints.
1603 include 'DIMENSIONS'
1604 include 'COMMON.IOUNITS'
1605 include 'COMMON.GEO'
1606 include 'COMMON.VAR'
1607 include 'COMMON.INTERACT'
1608 include 'COMMON.LOCAL'
1609 include 'COMMON.NAMES'
1610 include 'COMMON.CHAIN'
1611 include 'COMMON.FFIELD'
1612 include 'COMMON.SBRIDGE'
1613 include 'COMMON.HEADER'
1614 include 'COMMON.CONTROL'
1615 include 'COMMON.DBASE'
1616 include 'COMMON.THREAD'
1617 include 'COMMON.TIME1'
1618 integer i,j,itype_pdb(maxres)
1619 common /pizda/ itype_pdb
1620 double precision dist
1622 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1623 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1624 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1626 do i=nstart_sup,nstart_sup+nsup-1
1627 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1628 cd & ' seq_pdb', restyp(itype_pdb(i))
1629 do j=i+2,nstart_sup+nsup-1
1631 ihpb(nhpb)=i+nstart_seq-nstart_sup
1632 jhpb(nhpb)=j+nstart_seq-nstart_sup
1634 dhpb(nhpb)=dist(i,j)
1637 cd write (iout,'(a)') 'Distance constraints:'
1642 cd if (ii.gt.nres) then
1647 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1648 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1649 cd & dhpb(i),forcon(i)
1653 c----------------------------------------------------------------------------
1656 include 'DIMENSIONS'
1657 include 'COMMON.MAP'
1658 include 'COMMON.IOUNITS'
1660 character*3 angid(4) /'THE','PHI','ALP','OME'/
1661 character*80 mapcard,ucase
1663 read (inp,'(a)') mapcard
1664 mapcard=ucase(mapcard)
1665 if (index(mapcard,'PHI').gt.0) then
1667 else if (index(mapcard,'THE').gt.0) then
1669 else if (index(mapcard,'ALP').gt.0) then
1671 else if (index(mapcard,'OME').gt.0) then
1674 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1675 stop 'Error - illegal variable spec in MAP card.'
1677 call readi (mapcard,'RES1',res1(imap),0)
1678 call readi (mapcard,'RES2',res2(imap),0)
1679 if (res1(imap).eq.0) then
1680 res1(imap)=res2(imap)
1681 else if (res2(imap).eq.0) then
1682 res2(imap)=res1(imap)
1684 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1686 & 'Error - illegal definition of variable group in MAP.'
1687 stop 'Error - illegal definition of variable group in MAP.'
1689 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1690 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1691 call readi(mapcard,'NSTEP',nstep(imap),0)
1692 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1694 & 'Illegal boundary and/or step size specification in MAP.'
1695 stop 'Illegal boundary and/or step size specification in MAP.'
1700 c----------------------------------------------------------------------------
1703 include 'DIMENSIONS'
1704 include 'COMMON.IOUNITS'
1705 include 'COMMON.GEO'
1706 include 'COMMON.CSA'
1707 include 'COMMON.BANK'
1708 include 'COMMON.CONTROL'
1710 character*620 mcmcard
1711 call card_concat(mcmcard)
1713 call readi(mcmcard,'NCONF',nconf,50)
1714 call readi(mcmcard,'NADD',nadd,0)
1715 call readi(mcmcard,'JSTART',jstart,1)
1716 call readi(mcmcard,'JEND',jend,1)
1717 call readi(mcmcard,'NSTMAX',nstmax,500000)
1718 call readi(mcmcard,'N0',n0,1)
1719 call readi(mcmcard,'N1',n1,6)
1720 call readi(mcmcard,'N2',n2,4)
1721 call readi(mcmcard,'N3',n3,0)
1722 call readi(mcmcard,'N4',n4,0)
1723 call readi(mcmcard,'N5',n5,0)
1724 call readi(mcmcard,'N6',n6,10)
1725 call readi(mcmcard,'N7',n7,0)
1726 call readi(mcmcard,'N8',n8,0)
1727 call readi(mcmcard,'N9',n9,0)
1728 call readi(mcmcard,'N14',n14,0)
1729 call readi(mcmcard,'N15',n15,0)
1730 call readi(mcmcard,'N16',n16,0)
1731 call readi(mcmcard,'N17',n17,0)
1732 call readi(mcmcard,'N18',n18,0)
1734 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1736 call readi(mcmcard,'NDIFF',ndiff,2)
1737 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1738 call readi(mcmcard,'IS1',is1,1)
1739 call readi(mcmcard,'IS2',is2,8)
1740 call readi(mcmcard,'NRAN0',nran0,4)
1741 call readi(mcmcard,'NRAN1',nran1,2)
1742 call readi(mcmcard,'IRR',irr,1)
1743 call readi(mcmcard,'NSEED',nseed,20)
1744 call readi(mcmcard,'NTOTAL',ntotal,10000)
1745 call reada(mcmcard,'CUT1',cut1,2.0d0)
1746 call reada(mcmcard,'CUT2',cut2,5.0d0)
1747 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1748 call readi(mcmcard,'ICMAX',icmax,3)
1749 call readi(mcmcard,'IRESTART',irestart,0)
1750 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1753 call reada(mcmcard,'DELE',dele,20.0d0)
1754 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1755 call readi(mcmcard,'IREF',iref,0)
1756 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1757 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1758 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1759 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1760 write (iout,*) "NCONF_IN",nconf_in
1763 c----------------------------------------------------------------------------
1766 include 'DIMENSIONS'
1767 include 'COMMON.MCM'
1768 include 'COMMON.MCE'
1769 include 'COMMON.IOUNITS'
1771 character*320 mcmcard
1773 call card_concat(mcmcard)
1774 call readi(mcmcard,'MAXACC',maxacc,100)
1775 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1776 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1777 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1778 call readi(mcmcard,'MAXREPM',maxrepm,200)
1779 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1780 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1781 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1782 call reada(mcmcard,'E_UP',e_up,5.0D0)
1783 call reada(mcmcard,'DELTE',delte,0.1D0)
1784 call readi(mcmcard,'NSWEEP',nsweep,5)
1785 call readi(mcmcard,'NSTEPH',nsteph,0)
1786 call readi(mcmcard,'NSTEPC',nstepc,0)
1787 call reada(mcmcard,'TMIN',tmin,298.0D0)
1788 call reada(mcmcard,'TMAX',tmax,298.0D0)
1789 call readi(mcmcard,'NWINDOW',nwindow,0)
1790 call readi(mcmcard,'PRINT_MC',print_mc,0)
1791 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1792 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1793 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1794 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1795 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1796 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1797 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1798 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1799 if (nwindow.gt.0) then
1800 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1802 winlen(i)=winend(i)-winstart(i)+1
1805 if (tmax.lt.tmin) tmax=tmin
1806 if (tmax.eq.tmin) then
1810 if (nstepc.gt.0 .and. nsteph.gt.0) then
1811 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1812 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1814 C Probabilities of different move types
1815 sumpro_type(0)=0.0D0
1816 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1817 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1818 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1819 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1820 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1821 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1822 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1824 print *,'i',i,' sumprotype',sumpro_type(i)
1825 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1826 print *,'i',i,' sumprotype',sumpro_type(i)
1830 c----------------------------------------------------------------------------
1831 subroutine read_minim
1833 include 'DIMENSIONS'
1834 include 'COMMON.MINIM'
1835 include 'COMMON.IOUNITS'
1837 character*320 minimcard
1838 call card_concat(minimcard)
1839 call readi(minimcard,'MAXMIN',maxmin,2000)
1840 call readi(minimcard,'MAXFUN',maxfun,5000)
1841 call readi(minimcard,'MINMIN',minmin,maxmin)
1842 call readi(minimcard,'MINFUN',minfun,maxmin)
1843 call reada(minimcard,'TOLF',tolf,1.0D-2)
1844 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1845 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1846 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1847 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1848 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1849 & 'Options in energy minimization:'
1850 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1851 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1852 & 'MinMin:',MinMin,' MinFun:',MinFun,
1853 & ' TolF:',TolF,' RTolF:',RTolF
1856 c----------------------------------------------------------------------------
1857 subroutine read_angles(kanal,*)
1859 include 'DIMENSIONS'
1860 include 'COMMON.GEO'
1861 include 'COMMON.VAR'
1862 include 'COMMON.CHAIN'
1863 include 'COMMON.IOUNITS'
1864 include 'COMMON.CONTROL'
1866 c Read angles from input
1868 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1869 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1870 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1871 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1874 c 9/7/01 avoid 180 deg valence angle
1875 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1877 theta(i)=deg2rad*theta(i)
1878 phi(i)=deg2rad*phi(i)
1879 alph(i)=deg2rad*alph(i)
1880 omeg(i)=deg2rad*omeg(i)
1885 c----------------------------------------------------------------------------
1886 subroutine reada(rekord,lancuch,wartosc,default)
1888 character*(*) rekord,lancuch
1889 double precision wartosc,default
1892 iread=index(rekord,lancuch)
1893 if (iread.eq.0) then
1897 iread=iread+ilen(lancuch)+1
1898 read (rekord(iread:),*,err=10,end=10) wartosc
1903 c----------------------------------------------------------------------------
1904 subroutine readi(rekord,lancuch,wartosc,default)
1906 character*(*) rekord,lancuch
1907 integer wartosc,default
1910 iread=index(rekord,lancuch)
1911 if (iread.eq.0) then
1915 iread=iread+ilen(lancuch)+1
1916 read (rekord(iread:),*,err=10,end=10) wartosc
1921 c----------------------------------------------------------------------------
1922 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1925 integer tablica(dim),default
1926 character*(*) rekord,lancuch
1933 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1934 if (iread.eq.0) return
1935 iread=iread+ilen(lancuch)+1
1936 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1939 c----------------------------------------------------------------------------
1940 subroutine multreada(rekord,lancuch,tablica,dim,default)
1943 double precision tablica(dim),default
1944 character*(*) rekord,lancuch
1951 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1952 if (iread.eq.0) return
1953 iread=iread+ilen(lancuch)+1
1954 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1957 c----------------------------------------------------------------------------
1958 subroutine openunits
1960 include 'DIMENSIONS'
1964 character*16 form,nodename
1967 include 'COMMON.SETUP'
1968 include 'COMMON.IOUNITS'
1970 include 'COMMON.CONTROL'
1971 integer lenpre,lenpot,ilen,lentmp,npos
1973 character*3 out1file_text,ucase
1978 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1979 call getenv_loc("PREFIX",prefix)
1981 call getenv_loc("POT",pot)
1982 call getenv_loc("DIRTMP",tmpdir)
1983 call getenv_loc("CURDIR",curdir)
1984 call getenv_loc("OUT1FILE",out1file_text)
1985 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1986 out1file_text=ucase(out1file_text)
1987 if (out1file_text(1:1).eq."Y") then
1990 out1file=fg_rank.gt.0
1995 if (lentmp.gt.0) then
1996 write (*,'(80(1h!))')
1997 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1998 write (*,'(80(1h!))')
1999 write (*,*)"All output files will be on node /tmp directory."
2001 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2002 if (me.eq.king) then
2003 write (*,*) "The master node is ",nodename
2004 else if (fg_rank.eq.0) then
2005 write (*,*) "I am the CG slave node ",nodename
2007 write (*,*) "I am the FG slave node ",nodename
2010 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2011 lenpre = lentmp+lenpre+1
2013 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2014 C Get the names and open the input files
2015 #if defined(WINIFL) || defined(WINPGI)
2016 open(1,file=pref_orig(:ilen(pref_orig))//
2017 & '.inp',status='old',readonly,shared)
2018 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2019 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2020 C Get parameter filenames and open the parameter files.
2021 call getenv_loc('BONDPAR',bondname)
2022 open (ibond,file=bondname,status='old',readonly,shared)
2023 call getenv_loc('THETPAR',thetname)
2024 open (ithep,file=thetname,status='old',readonly,shared)
2025 call getenv_loc('ROTPAR',rotname)
2026 open (irotam,file=rotname,status='old',readonly,shared)
2027 call getenv_loc('TORPAR',torname)
2028 open (itorp,file=torname,status='old',readonly,shared)
2029 call getenv_loc('TORDPAR',tordname)
2030 open (itordp,file=tordname,status='old',readonly,shared)
2031 call getenv_loc('FOURIER',fouriername)
2032 open (ifourier,file=fouriername,status='old',readonly,shared)
2033 call getenv_loc('ELEPAR',elename)
2034 open (ielep,file=elename,status='old',readonly,shared)
2035 call getenv_loc('SIDEPAR',sidename)
2036 open (isidep,file=sidename,status='old',readonly,shared)
2037 call getenv_loc('LIPTRANPAR',liptranname)
2038 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2039 #elif (defined CRAY) || (defined AIX)
2040 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2042 c print *,"Processor",myrank," opened file 1"
2043 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2044 c print *,"Processor",myrank," opened file 9"
2045 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2046 C Get parameter filenames and open the parameter files.
2047 call getenv_loc('BONDPAR',bondname)
2048 open (ibond,file=bondname,status='old',action='read')
2049 c print *,"Processor",myrank," opened file IBOND"
2050 call getenv_loc('THETPAR',thetname)
2051 open (ithep,file=thetname,status='old',action='read')
2053 call getenv_loc('THETPARPDB',thetname_pdb)
2054 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2056 c print *,"Processor",myrank," opened file ITHEP"
2057 call getenv_loc('ROTPAR',rotname)
2058 open (irotam,file=rotname,status='old',action='read')
2060 call getenv_loc('ROTPARPDB',rotname_pdb)
2061 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2063 c print *,"Processor",myrank," opened file IROTAM"
2064 call getenv_loc('TORPAR',torname)
2065 open (itorp,file=torname,status='old',action='read')
2066 c print *,"Processor",myrank," opened file ITORP"
2067 call getenv_loc('TORDPAR',tordname)
2068 open (itordp,file=tordname,status='old',action='read')
2069 c print *,"Processor",myrank," opened file ITORDP"
2070 call getenv_loc('SCCORPAR',sccorname)
2071 open (isccor,file=sccorname,status='old',action='read')
2072 c print *,"Processor",myrank," opened file ISCCOR"
2073 call getenv_loc('FOURIER',fouriername)
2074 open (ifourier,file=fouriername,status='old',action='read')
2075 c print *,"Processor",myrank," opened file IFOURIER"
2076 call getenv_loc('ELEPAR',elename)
2077 open (ielep,file=elename,status='old',action='read')
2078 c print *,"Processor",myrank," opened file IELEP"
2079 call getenv_loc('SIDEPAR',sidename)
2080 open (isidep,file=sidename,status='old',action='read')
2081 call getenv_loc('LIPTRANPAR',liptranname)
2082 open (iliptranpar,file=liptranname,status='old',action='read')
2083 c print *,"Processor",myrank," opened file ISIDEP"
2084 c print *,"Processor",myrank," opened parameter files"
2086 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2087 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2088 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2089 C Get parameter filenames and open the parameter files.
2090 call getenv_loc('BONDPAR',bondname)
2091 open (ibond,file=bondname,status='old')
2092 call getenv_loc('THETPAR',thetname)
2093 open (ithep,file=thetname,status='old')
2095 call getenv_loc('THETPARPDB',thetname_pdb)
2096 open (ithep_pdb,file=thetname_pdb,status='old')
2098 call getenv_loc('ROTPAR',rotname)
2099 open (irotam,file=rotname,status='old')
2101 call getenv_loc('ROTPARPDB',rotname_pdb)
2102 open (irotam_pdb,file=rotname_pdb,status='old')
2104 call getenv_loc('TORPAR',torname)
2105 open (itorp,file=torname,status='old')
2106 call getenv_loc('TORDPAR',tordname)
2107 open (itordp,file=tordname,status='old')
2108 call getenv_loc('SCCORPAR',sccorname)
2109 open (isccor,file=sccorname,status='old')
2110 call getenv_loc('FOURIER',fouriername)
2111 open (ifourier,file=fouriername,status='old')
2112 call getenv_loc('ELEPAR',elename)
2113 open (ielep,file=elename,status='old')
2114 call getenv_loc('SIDEPAR',sidename)
2115 open (isidep,file=sidename,status='old')
2116 call getenv_loc('LIPTRANPAR',liptranname)
2117 open (iliptranpar,file=liptranname,status='old')
2119 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2121 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2122 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2123 C Get parameter filenames and open the parameter files.
2124 call getenv_loc('BONDPAR',bondname)
2125 open (ibond,file=bondname,status='old',readonly)
2126 call getenv_loc('THETPAR',thetname)
2127 open (ithep,file=thetname,status='old',readonly)
2128 call getenv_loc('ROTPAR',rotname)
2129 open (irotam,file=rotname,status='old',readonly)
2130 call getenv_loc('TORPAR',torname)
2131 open (itorp,file=torname,status='old',readonly)
2132 call getenv_loc('TORDPAR',tordname)
2133 open (itordp,file=tordname,status='old',readonly)
2134 call getenv_loc('SCCORPAR',sccorname)
2135 open (isccor,file=sccorname,status='old',readonly)
2137 call getenv_loc('THETPARPDB',thetname_pdb)
2138 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2140 call getenv_loc('FOURIER',fouriername)
2141 open (ifourier,file=fouriername,status='old',readonly)
2142 call getenv_loc('ELEPAR',elename)
2143 open (ielep,file=elename,status='old',readonly)
2144 call getenv_loc('SIDEPAR',sidename)
2145 open (isidep,file=sidename,status='old',readonly)
2146 call getenv_loc('LIPTRANPAR',liptranname)
2147 open (iliptranpar,file=liptranname,status='old',action='read')
2149 call getenv_loc('ROTPARPDB',rotname_pdb)
2150 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2154 call getenv_loc('TUBEPAR',tubename)
2155 #if defined(WINIFL) || defined(WINPGI)
2156 open (itube,file=tubename,status='old',readonly,shared)
2157 #elif (defined CRAY) || (defined AIX)
2158 open (itube,file=tubename,status='old',action='read')
2160 open (itube,file=tubename,status='old')
2162 open (itube,file=tubename,status='old',readonly)
2167 C 8/9/01 In the newest version SCp interaction constants are read from a file
2168 C Use -DOLDSCP to use hard-coded constants instead.
2170 call getenv_loc('SCPPAR',scpname)
2171 #if defined(WINIFL) || defined(WINPGI)
2172 open (iscpp,file=scpname,status='old',readonly,shared)
2173 #elif (defined CRAY) || (defined AIX)
2174 open (iscpp,file=scpname,status='old',action='read')
2176 open (iscpp,file=scpname,status='old')
2178 open (iscpp,file=scpname,status='old',readonly)
2181 call getenv_loc('PATTERN',patname)
2182 #if defined(WINIFL) || defined(WINPGI)
2183 open (icbase,file=patname,status='old',readonly,shared)
2184 #elif (defined CRAY) || (defined AIX)
2185 open (icbase,file=patname,status='old',action='read')
2187 open (icbase,file=patname,status='old')
2189 open (icbase,file=patname,status='old',readonly)
2192 C Open output file only for CG processes
2193 c print *,"Processor",myrank," fg_rank",fg_rank
2194 if (fg_rank.eq.0) then
2196 if (nodes.eq.1) then
2199 npos = dlog10(dfloat(nodes-1))+1
2201 if (npos.lt.3) npos=3
2202 write (liczba,'(i1)') npos
2203 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2205 write (liczba,form) me
2206 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2207 & liczba(:ilen(liczba))
2208 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2210 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2212 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2213 & liczba(:ilen(liczba))//'.mol2'
2214 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2215 & liczba(:ilen(liczba))//'.stat'
2217 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2218 & //liczba(:ilen(liczba))//'.stat')
2219 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2222 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2223 & liczba(:ilen(liczba))//'.const'
2228 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2229 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2230 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2231 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2232 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2234 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2236 rest2name=prefix(:ilen(prefix))//'.rst'
2238 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2241 #if defined(AIX) || defined(PGI) || defined(CRAY)
2242 if (me.eq.king .or. .not. out1file) then
2243 open(iout,file=outname,status='unknown')
2245 open(iout,file="/dev/null",status="unknown")
2249 if (fg_rank.gt.0) then
2250 write (liczba,'(i3.3)') myrank/nfgtasks
2251 write (ll,'(bz,i3.3)') fg_rank
2252 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2258 open(igeom,file=intname,status='unknown',position='append')
2259 open(ipdb,file=pdbname,status='unknown')
2260 open(imol2,file=mol2name,status='unknown')
2261 open(istat,file=statname,status='unknown',position='append')
2263 c1out open(iout,file=outname,status='unknown')
2266 if (me.eq.king .or. .not.out1file)
2267 & open(iout,file=outname,status='unknown')
2269 if (fg_rank.gt.0) then
2270 write (liczba,'(i3.3)') myrank/nfgtasks
2271 write (ll,'(bz,i3.3)') fg_rank
2272 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2277 open(igeom,file=intname,status='unknown',access='append')
2278 open(ipdb,file=pdbname,status='unknown')
2279 open(imol2,file=mol2name,status='unknown')
2280 open(istat,file=statname,status='unknown',access='append')
2282 c1out open(iout,file=outname,status='unknown')
2285 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2286 csa_seed=prefix(:lenpre)//'.CSA.seed'
2287 csa_history=prefix(:lenpre)//'.CSA.history'
2288 csa_bank=prefix(:lenpre)//'.CSA.bank'
2289 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2290 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2291 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2292 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2293 csa_int=prefix(:lenpre)//'.int'
2294 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2295 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2296 csa_in=prefix(:lenpre)//'.CSA.in'
2297 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2300 write (iout,'(80(1h-))')
2301 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2302 write (iout,'(80(1h-))')
2303 write (iout,*) "Input file : ",
2304 & pref_orig(:ilen(pref_orig))//'.inp'
2305 write (iout,*) "Output file : ",
2306 & outname(:ilen(outname))
2308 write (iout,*) "Sidechain potential file : ",
2309 & sidename(:ilen(sidename))
2311 write (iout,*) "SCp potential file : ",
2312 & scpname(:ilen(scpname))
2314 write (iout,*) "Electrostatic potential file : ",
2315 & elename(:ilen(elename))
2316 write (iout,*) "Cumulant coefficient file : ",
2317 & fouriername(:ilen(fouriername))
2318 write (iout,*) "Torsional parameter file : ",
2319 & torname(:ilen(torname))
2320 write (iout,*) "Double torsional parameter file : ",
2321 & tordname(:ilen(tordname))
2322 write (iout,*) "SCCOR parameter file : ",
2323 & sccorname(:ilen(sccorname))
2324 write (iout,*) "Bond & inertia constant file : ",
2325 & bondname(:ilen(bondname))
2326 write (iout,*) "Bending-torsion parameter file : ",
2327 & thetname(:ilen(thetname))
2328 write (iout,*) "Rotamer parameter file : ",
2329 & rotname(:ilen(rotname))
2330 write (iout,*) "Thetpdb parameter file : ",
2331 & thetname_pdb(:ilen(thetname_pdb))
2332 write (iout,*) "Threading database : ",
2333 & patname(:ilen(patname))
2335 &write (iout,*)" DIRTMP : ",
2337 write (iout,'(80(1h-))')
2341 c----------------------------------------------------------------------------
2342 subroutine card_concat(card)
2343 implicit real*8 (a-h,o-z)
2344 include 'DIMENSIONS'
2345 include 'COMMON.IOUNITS'
2347 character*80 karta,ucase
2349 read (inp,'(a)') karta
2352 do while (karta(80:80).eq.'&')
2353 card=card(:ilen(card)+1)//karta(:79)
2354 read (inp,'(a)') karta
2357 card=card(:ilen(card)+1)//karta
2360 c------------------------------------------------------------------------------
2363 include 'DIMENSIONS'
2364 include 'COMMON.CHAIN'
2365 include 'COMMON.IOUNITS'
2366 include 'COMMON.CONTROL'
2368 include 'COMMON.QRESTR'
2370 include 'COMMON.LAGRANGE.5diag'
2372 include 'COMMON.LAGRANGE'
2375 open(irest2,file=rest2name,status='unknown')
2376 read(irest2,*) totT,EK,potE,totE,t_bath
2379 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2382 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2385 read (irest2,*) iset
2390 c------------------------------------------------------------------------------
2391 subroutine read_fragments
2393 include 'DIMENSIONS'
2397 include 'COMMON.SETUP'
2398 include 'COMMON.CHAIN'
2399 include 'COMMON.IOUNITS'
2401 include 'COMMON.QRESTR'
2402 include 'COMMON.CONTROL'
2404 read(inp,*) nset,nfrag,npair,nfrag_back
2405 loc_qlike=(nfrag_back.lt.0)
2406 nfrag_back=iabs(nfrag_back)
2407 if(me.eq.king.or..not.out1file)
2408 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2409 & " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2411 read(inp,*) mset(iset)
2413 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2415 if(me.eq.king.or..not.out1file)
2416 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2417 & ifrag(2,i,iset), qinfrag(i,iset)
2420 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2422 if(me.eq.king.or..not.out1file)
2423 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2424 & ipair(2,i,iset), qinpair(i,iset)
2428 read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2429 & wfrag_back(2,i,iset),qin_back(2,i,iset),
2430 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2431 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2432 if(me.eq.king.or..not.out1file)
2433 & write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2434 & wfrag_back(2,i,iset),qin_back(3,i,iset),
2435 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2436 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2440 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2441 & wfrag_back(3,i,iset),
2442 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2443 if(me.eq.king.or..not.out1file)
2444 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2445 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2451 C---------------------------------------------------------------------------
2452 subroutine read_afminp
2454 include 'DIMENSIONS'
2458 include 'COMMON.SETUP'
2459 include 'COMMON.CONTROL'
2460 include 'COMMON.CHAIN'
2461 include 'COMMON.IOUNITS'
2462 include 'COMMON.SBRIDGE'
2463 character*320 afmcard
2465 c print *, "wchodze"
2466 call card_concat(afmcard)
2467 call readi(afmcard,"BEG",afmbeg,0)
2468 call readi(afmcard,"END",afmend,0)
2469 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2470 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2471 print *,'FORCE=' ,forceAFMconst
2472 CCCC NOW PROPERTIES FOR AFM
2475 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2477 distafminit=dsqrt(distafminit)
2478 c print *,'initdist',distafminit
2481 c-------------------------------------------------------------------------------
2482 subroutine read_saxs_constr
2484 include 'DIMENSIONS'
2488 include 'COMMON.SETUP'
2489 include 'COMMON.CONTROL'
2490 include 'COMMON.SAXS'
2491 include 'COMMON.CHAIN'
2492 include 'COMMON.IOUNITS'
2493 include 'COMMON.SBRIDGE'
2494 double precision cm(3),cnorm
2497 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2499 if (saxs_mode.eq.0) then
2500 c SAXS distance distribution
2502 read(inp,*) distsaxs(i),Psaxs(i)
2506 Cnorm = Cnorm + Psaxs(i)
2508 write (iout,*) "Cnorm",Cnorm
2510 Psaxs(i)=Psaxs(i)/Cnorm
2512 write (iout,*) "Normalized distance distribution from SAXS"
2514 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2518 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2520 write (iout,*) "Wsaxs0",Wsaxs0
2524 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2531 cm(j)=cm(j)+Csaxs(j,i)
2539 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2542 write (iout,*) "SAXS sphere coordinates"
2544 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2550 c-------------------------------------------------------------------------------
2551 subroutine read_dist_constr
2553 include 'DIMENSIONS'
2557 include 'COMMON.SETUP'
2558 include 'COMMON.CONTROL'
2559 include 'COMMON.CHAIN'
2560 include 'COMMON.IOUNITS'
2561 include 'COMMON.SBRIDGE'
2562 include 'COMMON.INTERACT'
2563 integer i,j,k,ii,jj,itemp,link_type,iiend,jjend,kk
2564 integer nfrag_,npair_,ndist_,ifrag_(2,100),ipair_(2,1000)
2565 double precision wfrag_(100),wpair_(1000)
2566 double precision ddjk,dist,dist_cut,fordepthmax
2567 character*5000 controlcard
2568 logical normalize,next
2570 double precision scal_bfac
2571 double precision xlink(4,0:4) /
2573 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
2574 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
2575 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
2576 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
2577 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
2578 c print *, "WCHODZE"
2579 write (iout,*) "Calling read_dist_constr"
2580 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2582 restr_on_coord=.false.
2591 call card_concat(controlcard)
2592 next = index(controlcard,"NEXT").gt.0
2593 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2594 write (iout,*) "restr_type",restr_type
2595 call readi(controlcard,"NFRAG",nfrag_,0)
2596 call readi(controlcard,"NFRAG",nfrag_,0)
2597 call readi(controlcard,"NPAIR",npair_,0)
2598 call readi(controlcard,"NDIST",ndist_,0)
2599 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2600 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2601 if (restr_type.eq.10)
2602 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2603 if (restr_type.eq.12)
2604 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2605 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2606 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2607 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2608 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2609 normalize = index(controlcard,"NORMALIZE").gt.0
2610 write (iout,*) "WBOLTZD",wboltzd
2611 write (iout,*) "SCAL_PEAK",scal_peak
2612 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2613 write (iout,*) "IFRAG"
2615 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2617 write (iout,*) "IPAIR"
2619 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2621 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2623 & "Distance restraints as generated from reference structure"
2625 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2626 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2627 & ifrag_(2,i)=nstart_sup+nsup-1
2628 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2630 if (wfrag_(i).eq.0.0d0) cycle
2631 do j=ifrag_(1,i),ifrag_(2,i)-1
2632 do k=j+1,ifrag_(2,i)
2633 c write (iout,*) "j",j," k",k
2635 if (restr_type.eq.1) then
2641 forcon(nhpb)=wfrag_(i)
2642 else if (constr_dist.eq.2) then
2643 if (ddjk.le.dist_cut) then
2649 forcon(nhpb)=wfrag_(i)
2651 else if (restr_type.eq.3) then
2657 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2660 if (.not.out1file .or. me.eq.king)
2661 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2662 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2664 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2665 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2671 if (wpair_(i).eq.0.0d0) cycle
2679 do j=ifrag_(1,ii),ifrag_(2,ii)
2680 do k=ifrag_(1,jj),ifrag_(2,jj)
2682 if (restr_type.eq.1) then
2688 forcon(nhpb)=wpair_(i)
2689 else if (constr_dist.eq.2) then
2690 if (ddjk.le.dist_cut) then
2696 forcon(nhpb)=wpair_(i)
2698 else if (restr_type.eq.3) then
2704 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2707 if (.not.out1file .or. me.eq.king)
2708 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2709 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2711 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2712 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2719 write (iout,*) "Distance restraints as read from input"
2721 if (restr_type.eq.12) then
2722 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2723 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2724 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2725 & fordepth_peak(nhpb_peak+1),npeak
2726 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2727 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2728 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2729 c & fordepth_peak(nhpb_peak+1),npeak
2730 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2731 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2732 nhpb_peak=nhpb_peak+1
2733 irestr_type_peak(nhpb_peak)=12
2734 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2737 if (.not.out1file .or. me.eq.king)
2738 & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2739 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2740 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2741 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2742 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2744 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2745 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2746 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2747 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2748 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2750 if (ibecarb_peak(nhpb_peak).eq.3) then
2751 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2752 else if (ibecarb_peak(nhpb_peak).eq.2) then
2753 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2754 else if (ibecarb_peak(nhpb_peak).eq.1) then
2755 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2756 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2758 else if (restr_type.eq.11) then
2759 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2760 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2761 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2762 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2764 irestr_type(nhpb)=11
2766 if (.not.out1file .or. me.eq.king)
2767 & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2768 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2769 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2771 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2772 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2773 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2775 c if (ibecarb(nhpb).gt.0) then
2776 c ihpb(nhpb)=ihpb(nhpb)+nres
2777 c jhpb(nhpb)=jhpb(nhpb)+nres
2779 if (ibecarb(nhpb).eq.3) then
2780 ihpb(nhpb)=ihpb(nhpb)+nres
2781 else if (ibecarb(nhpb).eq.2) then
2782 ihpb(nhpb)=ihpb(nhpb)+nres
2783 else if (ibecarb(nhpb).eq.1) then
2784 ihpb(nhpb)=ihpb(nhpb)+nres
2785 jhpb(nhpb)=jhpb(nhpb)+nres
2787 else if (restr_type.eq.10) then
2788 c Cross-lonk Markov-like potential
2789 call card_concat(controlcard)
2790 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2791 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2793 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2794 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2795 if (index(controlcard,"ZL").gt.0) then
2797 else if (index(controlcard,"ADH").gt.0) then
2799 else if (index(controlcard,"PDH").gt.0) then
2801 else if (index(controlcard,"DSS").gt.0) then
2806 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2807 & xlink(1,link_type))
2808 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2809 & xlink(2,link_type))
2810 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2811 & xlink(3,link_type))
2812 call reada(controlcard,"SIGMA",forcon(nhpb+1),
2813 & xlink(4,link_type))
2814 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2815 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2816 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2817 if (forcon(nhpb+1).le.0.0d0 .or.
2818 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2820 irestr_type(nhpb)=10
2821 if (ibecarb(nhpb).eq.3) then
2822 jhpb(nhpb)=jhpb(nhpb)+nres
2823 else if (ibecarb(nhpb).eq.2) then
2824 ihpb(nhpb)=ihpb(nhpb)+nres
2825 else if (ibecarb(nhpb).eq.1) then
2826 ihpb(nhpb)=ihpb(nhpb)+nres
2827 jhpb(nhpb)=jhpb(nhpb)+nres
2830 if (.not.out1file .or. me.eq.king)
2831 & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2832 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2833 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2836 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2837 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2838 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2843 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2844 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2845 if (forcon(nhpb+1).gt.0.0d0) then
2847 if (dhpb1(nhpb).eq.0.0d0) then
2852 if (ibecarb(nhpb).eq.3) then
2853 jhpb(nhpb)=jhpb(nhpb)+nres
2854 else if (ibecarb(nhpb).eq.2) then
2855 ihpb(nhpb)=ihpb(nhpb)+nres
2856 else if (ibecarb(nhpb).eq.1) then
2857 ihpb(nhpb)=ihpb(nhpb)+nres
2858 jhpb(nhpb)=jhpb(nhpb)+nres
2860 if (dhpb(nhpb).eq.0.0d0)
2861 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2864 if (.not.out1file .or. me.eq.king)
2865 & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2866 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2868 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2869 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2872 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2873 C if (forcon(nhpb+1).gt.0.0d0) then
2875 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2878 if (restr_type.eq.4) then
2879 write (iout,*) "The BFAC array"
2881 write (iout,'(i5,f10.5)') i,bfac(i)
2884 if (itype(i).eq.ntyp1) cycle
2886 if (itype(j).eq.ntyp1) cycle
2887 if (itype(i).eq.10) then
2892 if (itype(j).eq.10) then
2902 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2905 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2906 ihpb(nhpb)=i+nres*ii
2907 jhpb(nhpb)=j+nres*jj
2908 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2910 if (.not.out1file .or. me.eq.king) then
2911 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2912 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2913 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2917 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2918 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2919 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2929 if (restr_type.eq.5) then
2930 restr_on_coord=.true.
2932 if (itype(i).eq.ntyp1) cycle
2933 bfac(i)=(scal_bfac/bfac(i))**2
2942 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
2943 & fordepthmax=fordepth(i)
2946 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
2951 c-------------------------------------------------------------------------------
2952 subroutine read_constr_homology
2954 include 'DIMENSIONS'
2958 include 'COMMON.SETUP'
2959 include 'COMMON.CONTROL'
2960 include 'COMMON.HOMOLOGY'
2961 include 'COMMON.CHAIN'
2962 include 'COMMON.IOUNITS'
2964 include 'COMMON.QRESTR'
2965 include 'COMMON.GEO'
2966 include 'COMMON.INTERACT'
2967 include 'COMMON.NAMES'
2969 c For new homol impl
2971 include 'COMMON.VAR'
2974 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2976 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2977 c & sigma_odl_temp(maxres,maxres,max_template)
2979 character*24 model_ki_dist, model_ki_angle
2980 character*500 controlcard
2981 integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
2982 & ik,iistart,iishift
2987 c FP - Nov. 2014 Temporary specifications for new vars
2989 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
2991 double precision, dimension (max_template,maxres) :: rescore
2992 double precision, dimension (max_template,maxres) :: rescore2
2993 double precision, dimension (max_template,maxres) :: rescore3
2994 double precision distal
2995 character*24 pdbfile,tpl_k_rescore
2996 c -----------------------------------------------------------------
2997 c Reading multiple PDB ref structures and calculation of retraints
2998 c not using pre-computed ones stored in files model_ki_{dist,angle}
3000 c -----------------------------------------------------------------
3003 c Alternative: reading from input
3004 call card_concat(controlcard)
3005 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
3006 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
3007 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
3008 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
3009 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
3010 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
3011 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3012 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3013 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3014 if(.not.read2sigma.and.start_from_model) then
3015 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3016 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3017 start_from_model=.false.
3019 if(start_from_model .and. (me.eq.king .or. .not. out1file))
3020 & write(iout,*) 'START_FROM_MODELS is ON'
3021 if(start_from_model .and. rest) then
3022 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3023 write(iout,*) 'START_FROM_MODELS is OFF'
3024 write(iout,*) 'remove restart keyword from input'
3027 if (homol_nset.gt.1)then
3028 call card_concat(controlcard)
3029 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3030 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3031 write(iout,*) "iset homology_weight "
3033 write(iout,*) i,waga_homology(i)
3036 iset=mod(kolor,homol_nset)+1
3039 waga_homology(1)=1.0
3042 cd write (iout,*) "nnt",nnt," nct",nct
3049 c write(iout,*) 'nnt=',nnt,'nct=',nct
3052 do k=1,constr_homology
3065 if (read_homol_frag) then
3066 call read_klapaucjusz
3069 do k=1,constr_homology
3071 read(inp,'(a)') pdbfile
3072 if(me.eq.king .or. .not. out1file)
3073 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3074 & pdbfile(:ilen(pdbfile))
3075 open(ipdbin,file=pdbfile,status='old',err=33)
3077 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3078 & pdbfile(:ilen(pdbfile))
3081 c print *,'Begin reading pdb data'
3083 c Files containing res sim or local scores (former containing sigmas)
3086 write(kic2,'(bz,i2.2)') k
3088 tpl_k_rescore="template"//kic2//".sco"
3091 if (read2sigma) then
3092 call readpdb_template(k)
3097 c Distance restraints
3100 C Copy the coordinates from reference coordinates (?)
3104 c write (iout,*) "c(",j,i,") =",c(j,i)
3108 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3110 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3111 open (ientin,file=tpl_k_rescore,status='old')
3112 if (nnt.gt.1) rescore(k,1)=0.0d0
3113 do irec=nnt,nct ! loop for reading res sim
3114 if (read2sigma) then
3115 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3116 & rescore3_tmp,idomain_tmp
3118 idomain(k,i_tmp)=idomain_tmp
3119 rescore(k,i_tmp)=rescore_tmp
3120 rescore2(k,i_tmp)=rescore2_tmp
3121 rescore3(k,i_tmp)=rescore3_tmp
3122 if (.not. out1file .or. me.eq.king)
3123 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3124 & i_tmp,rescore2_tmp,rescore_tmp,
3125 & rescore3_tmp,idomain_tmp
3128 read (ientin,*,end=1401) rescore_tmp
3130 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3131 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3132 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3137 if (waga_dist.ne.0.0d0) then
3145 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3146 c write (iout,*) k,i,j,distal,dist2_cut
3148 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3149 & .and. distal.le.dist2_cut ) then
3155 c write (iout,*) "k",k
3156 c write (iout,*) "i",i," j",j," constr_homology",
3161 if (read2sigma) then
3164 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3166 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3167 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3168 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3170 if (odl(k,ii).le.dist_cut) then
3171 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3174 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3175 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3177 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3178 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3182 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3185 l_homo(k,ii)=.false.
3192 c Theta, dihedral and SC retraints
3194 if (waga_angle.gt.0.0d0) then
3195 c open (ientin,file=tpl_k_sigma_dih,status='old')
3196 c do irec=1,maxres-3 ! loop for reading sigma_dih
3197 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3198 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3199 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3200 c & sigma_dih(k,i+nnt-1)
3205 if (idomain(k,i).eq.0) then
3209 dih(k,i)=phiref(i) ! right?
3210 c read (ientin,*) sigma_dih(k,i) ! original variant
3211 c write (iout,*) "dih(",k,i,") =",dih(k,i)
3212 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3213 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3214 c & "rescore(",k,i-2,") =",rescore(k,i-2),
3215 c & "rescore(",k,i-3,") =",rescore(k,i-3)
3217 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3218 & rescore(k,i-2)+rescore(k,i-3))/4.0
3219 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3220 c write (iout,*) "Raw sigmas for dihedral angle restraints"
3221 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3222 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3223 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
3224 c Instead of res sim other local measure of b/b str reliability possible
3225 if (sigma_dih(k,i).ne.0)
3226 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3227 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3232 if (waga_theta.gt.0.0d0) then
3233 c open (ientin,file=tpl_k_sigma_theta,status='old')
3234 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3235 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3236 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3237 c & sigma_theta(k,i+nnt-1)
3242 do i = nnt+2,nct ! right? without parallel.
3243 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3244 c do i=ithet_start,ithet_end ! with FG parallel.
3245 if (idomain(k,i).eq.0) then
3246 sigma_theta(k,i)=0.0
3249 thetatpl(k,i)=thetaref(i)
3250 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3251 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3252 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3253 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3254 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3255 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3256 & rescore(k,i-2))/3.0
3257 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3258 if (sigma_theta(k,i).ne.0)
3259 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3261 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3262 c rescore(k,i-2) ! right expression ?
3263 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3267 if (waga_d.gt.0.0d0) then
3268 c open (ientin,file=tpl_k_sigma_d,status='old')
3269 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3270 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3271 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3272 c & sigma_d(k,i+nnt-1)
3276 do i = nnt,nct ! right? without parallel.
3277 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3278 c do i=loc_start,loc_end ! with FG parallel.
3279 if (itype(i).eq.10) cycle
3280 if (idomain(k,i).eq.0 ) then
3287 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3288 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3289 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3290 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3291 sigma_d(k,i)=rescore3(k,i) ! right expression ?
3292 if (sigma_d(k,i).ne.0)
3293 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3295 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3296 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3297 c read (ientin,*) sigma_d(k,i) ! 1st variant
3302 c remove distance restraints not used in any model from the list
3303 c shift data in all arrays
3305 if (waga_dist.ne.0.0d0) then
3311 if (ii_in_use(ii).eq.0.and.liiflag) then
3315 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3316 & .not.liiflag.and.ii.eq.lim_odl) then
3317 if (ii.eq.lim_odl) then
3318 iishift=ii-iistart+1
3323 do ki=iistart,lim_odl-iishift
3324 ires_homo(ki)=ires_homo(ki+iishift)
3325 jres_homo(ki)=jres_homo(ki+iishift)
3326 ii_in_use(ki)=ii_in_use(ki+iishift)
3327 do k=1,constr_homology
3328 odl(k,ki)=odl(k,ki+iishift)
3329 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3330 l_homo(k,ki)=l_homo(k,ki+iishift)
3334 lim_odl=lim_odl-iishift
3340 endif ! .not. klapaucjusz
3342 if (constr_homology.gt.0) call homology_partition
3343 if (constr_homology.gt.0) call init_int_table
3344 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3345 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3349 if (.not.out_template_restr) return
3350 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3351 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3352 write (iout,*) "Distance restraints from templates"
3354 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3355 & ii,ires_homo(ii),jres_homo(ii),
3356 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3357 & ki=1,constr_homology)
3359 write (iout,*) "Dihedral angle restraints from templates"
3361 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3362 & (rad2deg*dih(ki,i),
3363 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3365 write (iout,*) "Virtual-bond angle restraints from templates"
3367 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3368 & (rad2deg*thetatpl(ki,i),
3369 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3371 write (iout,*) "SC restraints from templates"
3373 write(iout,'(i5,100(4f8.2,4x))') i,
3374 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3375 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3378 c -----------------------------------------------------------------
3381 c----------------------------------------------------------------------
3383 subroutine flush(iu)
3388 subroutine flush(iu)
3393 c------------------------------------------------------------------------------
3394 subroutine copy_to_tmp(source)
3396 include "DIMENSIONS"
3397 include "COMMON.IOUNITS"
3398 character*(*) source
3399 character* 256 tmpfile
3403 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3404 inquire(file=tmpfile,exist=ex)
3406 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3407 & " to temporary directory..."
3408 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3409 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3413 c------------------------------------------------------------------------------
3414 subroutine move_from_tmp(source)
3416 include "DIMENSIONS"
3417 include "COMMON.IOUNITS"
3418 character*(*) source
3421 write (*,*) "Moving ",source(:ilen(source)),
3422 & " from temporary directory to working directory"
3423 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3424 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3427 c------------------------------------------------------------------------------
3428 subroutine random_init(seed)
3430 C Initialize random number generator
3433 include 'DIMENSIONS'
3436 logical OKRandom, prng_restart
3438 integer iseed_array(4)
3439 integer error_msg,ierr
3441 include 'COMMON.IOUNITS'
3442 include 'COMMON.TIME1'
3443 include 'COMMON.THREAD'
3444 include 'COMMON.SBRIDGE'
3445 include 'COMMON.CONTROL'
3446 include 'COMMON.MCM'
3447 include 'COMMON.MAP'
3448 include 'COMMON.HEADER'
3449 include 'COMMON.CSA'
3450 include 'COMMON.CHAIN'
3451 include 'COMMON.MUCA'
3453 include 'COMMON.FFIELD'
3454 include 'COMMON.SETUP'
3456 double precision seed,ran_number
3457 iseed=-dint(dabs(seed))
3458 if (iseed.eq.0) then
3459 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3460 & 'Random seed undefined. The program will stop.'
3461 write (*,'(/80(1h*)/20x,a/80(1h*))')
3462 & 'Random seed undefined. The program will stop.'
3464 call mpi_finalize(mpi_comm_world,ierr)
3466 stop 'Bad random seed.'
3469 if (fg_rank.eq.0) then
3473 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3474 OKRandom = prng_restart(me,iseed)
3477 tmp=65536.0d0**(4-i)
3478 iseed_array(i) = dint(seed/tmp)
3479 seed=seed-iseed_array(i)*tmp
3482 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3483 & (iseed_array(i),i=1,4)
3484 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3485 & (iseed_array(i),i=1,4)
3486 OKRandom = prng_restart(me,iseed_array)
3489 c r1 = prng_next(me)
3490 r1=ran_number(0.0D0,1.0D0)
3492 & write (iout,*) 'ran_num',r1
3493 if (r1.lt.0.0d0) OKRandom=.false.
3495 if (.not.OKRandom) then
3496 write (iout,*) 'PRNG IS NOT WORKING!!!'
3497 print *,'PRNG IS NOT WORKING!!!'
3500 call mpi_abort(mpi_comm_world,error_msg,ierr)
3503 write (iout,*) 'too many processors for parallel prng'
3504 write (*,*) 'too many processors for parallel prng'
3512 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3516 c----------------------------------------------------------------------
3517 subroutine read_klapaucjusz
3519 include 'DIMENSIONS'
3523 include 'COMMON.SETUP'
3524 include 'COMMON.CONTROL'
3525 include 'COMMON.HOMOLOGY'
3526 include 'COMMON.CHAIN'
3527 include 'COMMON.IOUNITS'
3529 include 'COMMON.GEO'
3530 include 'COMMON.INTERACT'
3531 include 'COMMON.NAMES'
3532 character*256 fragfile
3533 integer ninclust(maxclust),inclust(max_template,maxclust),
3534 & nresclust(maxclust),iresclust(maxres,maxclust),nclust
3537 character*24 model_ki_dist, model_ki_angle
3538 character*500 controlcard
3539 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp,
3540 & ik,ll,ii,kk,iistart,iishift,lim_xx
3541 double precision distal
3542 logical lprn /.true./
3548 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3549 double precision, dimension (max_template,maxres) :: rescore
3550 double precision, dimension (max_template,maxres) :: rescore2
3551 character*24 pdbfile,tpl_k_rescore
3554 c For new homol impl
3556 include 'COMMON.VAR'
3558 call getenv("FRAGFILE",fragfile)
3559 open(ientin,file=fragfile,status="old",err=10)
3560 read(ientin,*) constr_homology,nclust
3566 do k=1,constr_homology
3567 read(ientin,'(a)') pdbfile
3568 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3569 & pdbfile(:ilen(pdbfile))
3570 open(ipdbin,file=pdbfile,status='old',err=33)
3572 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3573 & pdbfile(:ilen(pdbfile))
3577 call readpdb_template(k)
3585 read(ientin,*) ninclust(i),nresclust(i)
3586 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3587 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3590 c Loop over clusters
3593 do ll = 1,ninclust(l)
3601 idomain(k,iresclust(i,l)+1) = 1
3603 idomain(k,iresclust(i,l)) = 1
3607 c Distance restraints
3610 C Copy the coordinates from reference coordinates (?)
3614 c write (iout,*) "c(",j,i,") =",c(j,i)
3617 call int_from_cart(.true.,.false.)
3618 call sc_loc_geom(.false.)
3620 thetaref(i)=theta(i)
3623 if (waga_dist.ne.0.0d0) then
3631 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3632 c write (iout,*) k,i,j,distal,dist2_cut
3634 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3635 & .and. distal.le.dist2_cut ) then
3641 c write (iout,*) "k",k
3642 c write (iout,*) "i",i," j",j," constr_homology",
3647 if (read2sigma) then
3650 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3652 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3653 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3654 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3656 if (odl(k,ii).le.dist_cut) then
3657 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3660 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3661 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3663 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3664 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3668 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3671 c l_homo(k,ii)=.false.
3678 c Theta, dihedral and SC retraints
3680 if (waga_angle.gt.0.0d0) then
3682 if (idomain(k,i).eq.0) then
3683 c sigma_dih(k,i)=0.0
3687 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3688 & rescore(k,i-2)+rescore(k,i-3))/4.0
3689 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3690 c & " sigma_dihed",sigma_dih(k,i)
3691 if (sigma_dih(k,i).ne.0)
3692 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3697 if (waga_theta.gt.0.0d0) then
3699 if (idomain(k,i).eq.0) then
3700 c sigma_theta(k,i)=0.0
3703 thetatpl(k,i)=thetaref(i)
3704 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3705 & rescore(k,i-2))/3.0
3706 if (sigma_theta(k,i).ne.0)
3707 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3711 if (waga_d.gt.0.0d0) then
3713 if (itype(i).eq.10) cycle
3714 if (idomain(k,i).eq.0 ) then
3721 sigma_d(k,i)=rescore(k,i)
3722 if (sigma_d(k,i).ne.0)
3723 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3724 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3730 c remove distance restraints not used in any model from the list
3731 c shift data in all arrays
3733 if (waga_dist.ne.0.0d0) then
3739 if (ii_in_use(ii).eq.0.and.liiflag) then
3743 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3744 & .not.liiflag.and.ii.eq.lim_odl) then
3745 if (ii.eq.lim_odl) then
3746 iishift=ii-iistart+1
3751 do ki=iistart,lim_odl-iishift
3752 ires_homo(ki)=ires_homo(ki+iishift)
3753 jres_homo(ki)=jres_homo(ki+iishift)
3754 ii_in_use(ki)=ii_in_use(ki+iishift)
3755 do k=1,constr_homology
3756 odl(k,ki)=odl(k,ki+iishift)
3757 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3758 l_homo(k,ki)=l_homo(k,ki+iishift)
3762 lim_odl=lim_odl-iishift
3769 10 stop "Error in fragment file"