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,'SEARCHSC').gt.0 .and.
181 & index(controlcard,'NOSEARCHSC').eq.0)
182 sideadd=(index(controlcard,'SIDEADD').gt.0)
183 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
184 mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
185 outpdb=(index(controlcard,'PDBOUT').gt.0)
186 outmol2=(index(controlcard,'MOL2OUT').gt.0)
187 pdbref=(index(controlcard,'PDBREF').gt.0)
188 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
189 indpdb=index(controlcard,'PDBSTART')
190 AFMlog=(index(controlcard,'AFM'))
191 selfguide=(index(controlcard,'SELFGUIDE'))
192 c print *,'AFMlog',AFMlog,selfguide,"KUPA"
193 call readi(controlcard,'TUBEMOD',tubelog,0)
194 c write (iout,*) TUBElog,"TUBEMODE"
195 call readi(controlcard,'IPRINT',iprint,0)
196 C SHIELD keyword sets if the shielding effect of side-chains is used
197 C 0 denots no shielding is used all peptide are equally despite the
198 C solvent accesible area
199 C 1 the newly introduced function
200 C 2 reseved for further possible developement
201 call readi(controlcard,'SHIELD',shield_mode,0)
202 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
203 write(iout,*) "shield_mode",shield_mode
205 call readi(controlcard,'TORMODE',tor_mode,0)
206 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
207 write(iout,*) "torsional and valence angle mode",tor_mode
208 call readi(controlcard,'MAXGEN',maxgen,10000)
209 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
210 call readi(controlcard,"KDIAG",kdiag,0)
211 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
212 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
213 & write (iout,*) "RESCALE_MODE",rescale_mode
214 split_ene=index(controlcard,'SPLIT_ENE').gt.0
215 if (index(controlcard,'REGULAR').gt.0.0D0) then
216 call reada(controlcard,'WEIDIS',weidis,0.1D0)
221 if (index(controlcard,"CASC").gt.0) then
223 c else if (index(controlcard,"CAONLY").gt.0) then
225 else if (index(controlcard,"SCONLY").gt.0) then
229 c write (iout,*) "ERROR! Must specify CASC CAONLY or SCONLY when",
230 c & " specifying REFSTR, PDBREF or REGULAR."
234 if (index(controlcard,'CHECKGRAD').gt.0) then
236 if (index(controlcard,' CART').gt.0) then
238 elseif (index(controlcard,'CARINT').gt.0) then
243 call reada(controlcard,'DELTA',aincr,1.0d-7)
244 c write (iout,*) "icheckgrad",icheckgrad
245 elseif (index(controlcard,'THREAD').gt.0) then
247 call readi(controlcard,'THREAD',nthread,0)
248 if (nthread.gt.0) then
249 call reada(controlcard,'WEIDIS',weidis,0.1D0)
252 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
253 stop 'Error termination in Read_Control.'
255 else if (index(controlcard,'MCMA').gt.0) then
257 else if (index(controlcard,'MCEE').gt.0) then
259 else if (index(controlcard,'MULTCONF').gt.0) then
261 else if (index(controlcard,'MAP').gt.0) then
263 call readi(controlcard,'MAP',nmap,0)
264 else if (index(controlcard,'CSA').gt.0) then
266 crc else if (index(controlcard,'ZSCORE').gt.0) then
268 crc ZSCORE is rm from UNRES, modecalc=9 is available
271 cfcm else if (index(controlcard,'MCMF').gt.0) then
273 else if (index(controlcard,'SOFTREG').gt.0) then
275 else if (index(controlcard,'CHECK_BOND').gt.0) then
277 else if (index(controlcard,'TEST').gt.0) then
279 else if (index(controlcard,'MD').gt.0) then
281 else if (index(controlcard,'RE ').gt.0) then
285 lmuca=index(controlcard,'MUCA').gt.0
286 call readi(controlcard,'MUCADYN',mucadyn,0)
287 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
288 if (lmuca .and. (me.eq.king .or. .not.out1file ))
290 write (iout,*) 'MUCADYN=',mucadyn
291 write (iout,*) 'MUCASMOOTH=',muca_smooth
294 iscode=index(controlcard,'ONE_LETTER')
295 indphi=index(controlcard,'PHI')
296 indback=index(controlcard,'BACK')
297 iranconf=index(controlcard,'RAND_CONF')
298 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
299 extconf=(index(controlcard,'EXTCONF').gt.0)
300 if (start_from_model) then
304 i2ndstr=index(controlcard,'USE_SEC_PRED')
305 gradout=index(controlcard,'GRADOUT').gt.0
306 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
307 C DISTCHAINMAX become obsolete for periodic boundry condition
308 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
309 C Reading the dimensions of box in x,y,z coordinates
310 call reada(controlcard,'BOXX',boxxsize,100.0d0)
311 call reada(controlcard,'BOXY',boxysize,100.0d0)
312 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
313 write(iout,*) "Periodic box dimensions",boxxsize,boxysize,boxzsize
314 c Cutoff range for interactions
315 call reada(controlcard,"R_CUT_INT",r_cut_int,25.0d0)
316 call reada(controlcard,"R_CUT_RESPA",r_cut_respa,2.0d0)
317 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
318 write (iout,*) "Cutoff on interactions",r_cut_int
320 & "Cutoff in switching short and long range interactions in RESPA",
322 write (iout,*) "lambda in switch function",rlamb
323 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
324 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
325 if (lipthick.gt.0.0d0) then
326 bordliptop=(boxzsize+lipthick)/2.0
327 bordlipbot=bordliptop-lipthick
329 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
330 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
331 buflipbot=bordlipbot+lipbufthick
332 bufliptop=bordliptop-lipbufthick
333 if ((lipbufthick*2.0d0).gt.lipthick)
334 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
336 write(iout,*) "bordliptop=",bordliptop
337 write(iout,*) "bordlipbot=",bordlipbot
338 write(iout,*) "bufliptop=",bufliptop
339 write(iout,*) "buflipbot=",buflipbot
340 write (iout,*) "SHIELD MODE",shield_mode
341 if (TUBElog.gt.0) then
342 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
343 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
344 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
345 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
346 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
347 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
348 buftubebot=bordtubebot+tubebufthick
349 buftubetop=bordtubetop-tubebufthick
351 if (me.eq.king .or. .not.out1file )
352 & write (iout,*) "DISTCHAINMAX",distchainmax
354 if(me.eq.king.or..not.out1file)
355 & write (iout,'(2a)') diagmeth(kdiag),
356 & ' routine used to diagonalize matrices.'
359 c--------------------------------------------------------------------------
360 subroutine read_REMDpar
366 include 'COMMON.IOUNITS'
367 include 'COMMON.TIME1'
370 include 'COMMON.LANGEVIN'
373 include 'COMMON.LANGEVIN.lang0.5diag'
375 include 'COMMON.LANGEVIN.lang0'
378 include 'COMMON.INTERACT'
379 include 'COMMON.NAMES'
381 include 'COMMON.REMD'
382 include 'COMMON.CONTROL'
383 include 'COMMON.SETUP'
385 character*320 controlcard
386 character*3200 controlcard1
387 integer iremd_m_total,i
389 if(me.eq.king.or..not.out1file)
390 & write (iout,*) "REMD setup"
392 call card_concat(controlcard)
393 call readi(controlcard,"NREP",nrep,3)
394 call readi(controlcard,"NSTEX",nstex,1000)
395 call reada(controlcard,"RETMIN",retmin,10.0d0)
396 call reada(controlcard,"RETMAX",retmax,1000.0d0)
397 mremdsync=(index(controlcard,'SYNC').gt.0)
398 call readi(controlcard,"NSYN",i_sync_step,100)
399 restart1file=(index(controlcard,'REST1FILE').gt.0)
400 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
401 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
402 if(max_cache_traj_use.gt.max_cache_traj)
403 & max_cache_traj_use=max_cache_traj
404 if(me.eq.king.or..not.out1file) then
405 cd if (traj1file) then
406 crc caching is in testing - NTWX is not ignored
407 cd write (iout,*) "NTWX value is ignored"
408 cd write (iout,*) " trajectory is stored to one file by master"
409 cd write (iout,*) " before exchange at NSTEX intervals"
411 write (iout,*) "NREP= ",nrep
412 write (iout,*) "NSTEX= ",nstex
413 write (iout,*) "SYNC= ",mremdsync
414 write (iout,*) "NSYN= ",i_sync_step
415 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
418 if (index(controlcard,'TLIST').gt.0) then
420 call card_concat(controlcard1)
421 read(controlcard1,*) (remd_t(i),i=1,nrep)
422 if(me.eq.king.or..not.out1file)
423 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
426 if (index(controlcard,'MLIST').gt.0) then
428 call card_concat(controlcard1)
429 read(controlcard1,*) (remd_m(i),i=1,nrep)
430 if(me.eq.king.or..not.out1file) then
431 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
434 iremd_m_total=iremd_m_total+remd_m(i)
436 write (iout,*) 'Total number of replicas ',iremd_m_total
441 & "Adaptive (PMF-biased) umbrella sampling will be run"
444 if(me.eq.king.or..not.out1file)
445 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
448 c--------------------------------------------------------------------------
449 subroutine read_MDpar
455 include 'COMMON.IOUNITS'
456 include 'COMMON.TIME1'
458 include 'COMMON.QRESTR'
460 include 'COMMON.LANGEVIN'
463 include 'COMMON.LANGEVIN.lang0.5diag'
465 include 'COMMON.LANGEVIN.lang0'
468 include 'COMMON.INTERACT'
469 include 'COMMON.NAMES'
471 include 'COMMON.SETUP'
472 include 'COMMON.CONTROL'
473 include 'COMMON.SPLITELE'
474 include 'COMMON.FFIELD'
476 character*320 controlcard
480 call card_concat(controlcard)
481 call readi(controlcard,"NSTEP",n_timestep,1000000)
482 call readi(controlcard,"NTWE",ntwe,100)
483 call readi(controlcard,"NTWX",ntwx,1000)
484 call reada(controlcard,"DT",d_time,1.0d-1)
485 call reada(controlcard,"DVMAX",dvmax,2.0d1)
486 call reada(controlcard,"DAMAX",damax,1.0d1)
487 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
488 call readi(controlcard,"LANG",lang,0)
489 RESPA = index(controlcard,"RESPA") .gt. 0
490 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
491 ntime_split0=ntime_split
492 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
493 ntime_split0=ntime_split
494 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
495 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
496 rest = index(controlcard,"REST").gt.0
497 tbf = index(controlcard,"TBF").gt.0
498 usampl = index(controlcard,"USAMPL").gt.0
499 scale_umb = index(controlcard,"SCALE_UMB").gt.0
500 adaptive = index(controlcard,"ADAPTIVE").gt.0
501 mdpdb = index(controlcard,"MDPDB").gt.0
502 call reada(controlcard,"T_BATH",t_bath,300.0d0)
503 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
504 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
505 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
506 if (count_reset_moment.eq.0) count_reset_moment=1000000000
507 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
508 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
509 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
510 if (count_reset_vel.eq.0) count_reset_vel=1000000000
511 large = index(controlcard,"LARGE").gt.0
512 print_compon = index(controlcard,"PRINT_COMPON").gt.0
513 rattle = index(controlcard,"RATTLE").gt.0
514 preminim = index(controlcard,"PREMINIM").gt.0
515 if (iranconf.gt.0 .or. indpdb.gt.0 .or. start_from_model) then
519 write (iout,*) "PREMINIM ",preminim
521 if (index(controlcard,'CART').gt.0) dccart=.true.
522 write (iout,*) "dccart ",dccart
523 write (iout,*) "read_minim ",index(controlcard,'READ_MINIM_PAR')
524 if (index(controlcard,'READ_MINIM_PAR').gt.0) call read_minim
526 c if performing umbrella sampling, fragments constrained are read from the fragment file
529 write (iout,*) "Umbrella sampling will be run"
530 if (scale_umb.and.adaptive) then
531 write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
532 write (iout,*) "Select one of those and re-run the job."
535 if (scale_umb) write (iout,*)
536 &"Umbrella-restraint force constants will be scaled by temperature"
540 if(me.eq.king.or..not.out1file) then
542 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
544 write (iout,'(a)') "The units are:"
545 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
546 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
547 & " acceleration: angstrom/(48.9 fs)**2"
548 write (iout,'(a)') "energy: kcal/mol, temperature: K"
550 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
551 write (iout,'(a60,f10.5,a)')
552 & "Initial time step of numerical integration:",d_time,
554 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
555 write (iout,'(a60,f10.5,a)') "Cutoff on interactions",r_cut_int,
557 write(iout,'(a60,i5)')"Frequency of updating interaction list",
560 write (iout,'(2a,i4,a)')
561 & "A-MTS algorithm used; initial time step for fast-varying",
562 & " short-range forces split into",ntime_split," steps."
563 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
564 & r_cut_respa," lambda",rlamb
566 write (iout,'(2a,f10.5)')
567 & "Maximum acceleration threshold to reduce the time step",
568 & "/increase split number:",damax
569 write (iout,'(2a,f10.5)')
570 & "Maximum predicted energy drift to reduce the timestep",
571 & "/increase split number:",edriftmax
572 write (iout,'(a60,f10.5)')
573 & "Maximum velocity threshold to reduce velocities:",dvmax
574 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
575 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
576 if (rattle) write (iout,'(a60)')
577 & "Rattle algorithm used to constrain the virtual bonds"
581 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
582 call reada(controlcard,"RWAT",rwat,1.4d0)
583 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
584 surfarea=index(controlcard,"SURFAREA").gt.0
585 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
586 if(me.eq.king.or..not.out1file)then
587 write (iout,'(/a,$)') "Langevin dynamics calculation"
590 & " with direct integration of Langevin equations"
591 else if (lang.eq.2) then
592 write (iout,'(a/)') " with TINKER stochasic MD integrator"
593 else if (lang.eq.3) then
594 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
595 else if (lang.eq.4) then
596 write (iout,'(a/)') " in overdamped mode"
598 write (iout,'(//a,i5)')
599 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
602 write (iout,'(a60,f10.5)') "Temperature:",t_bath
603 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
604 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
605 write (iout,'(a60,f10.5)')
606 & "Scaling factor of the friction forces:",scal_fric
607 if (surfarea) write (iout,'(2a,i10,a)')
608 & "Friction coefficients will be scaled by solvent-accessible",
609 & " surface area every",reset_fricmat," steps."
611 c Calculate friction coefficients and bounds of stochastic forces
612 eta=6*pi*cPoise*etawat
613 if(me.eq.king.or..not.out1file)
614 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
616 gamp=scal_fric*(pstok+rwat)*eta
617 stdfp=dsqrt(2*Rb*t_bath/d_time)
619 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
620 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
622 if(me.eq.king.or..not.out1file)then
623 write (iout,'(/2a/)')
624 & "Radii of site types and friction coefficients and std's of",
625 & " stochastic forces of fully exposed sites"
626 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
628 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
629 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
633 if(me.eq.king.or..not.out1file)then
634 write (iout,'(a)') "Berendsen bath calculation"
635 write (iout,'(a60,f10.5)') "Temperature:",t_bath
636 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
638 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
639 & count_reset_moment," steps"
641 & write (iout,'(a,i10,a)')
642 & "Velocities will be reset at random every",count_reset_vel,
646 if(me.eq.king.or..not.out1file)
647 & write (iout,'(a31)') "Microcanonical mode calculation"
649 if(me.eq.king.or..not.out1file)then
650 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
652 write(iout,*) "MD running with constraints."
653 write(iout,*) "Equilibration time ", eq_time, " mtus."
654 write(iout,*) "Constraining ", nfrag," fragments."
655 write(iout,*) "Length of each fragment, weight and q0:"
657 write (iout,*) "Set of restraints #",iset
659 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
660 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
662 write(iout,*) "constraints between ", npair, "fragments."
663 write(iout,*) "constraint pairs, weights and q0:"
665 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
666 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
668 write(iout,*) "angle constraints within ", nfrag_back,
669 & "backbone fragments."
671 write(iout,*) "fragment, weights, q0:"
673 write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
674 & ifrag_back(2,i,iset),
675 & wfrag_back(1,i,iset),qin_back(1,i,iset),
676 & wfrag_back(2,i,iset),qin_back(2,i,iset),
677 & wfrag_back(3,i,iset),qin_back(3,i,iset)
680 write(iout,*) "fragment, weights:"
682 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
683 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
684 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
688 iset=mod(kolor,nset)+1
691 if(me.eq.king.or..not.out1file)
692 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
695 c------------------------------------------------------------------------------
698 C Read molecular data.
704 integer error_msg,ierror,ierr,ierrcode
706 include 'COMMON.IOUNITS'
709 include 'COMMON.INTERACT'
710 include 'COMMON.LOCAL'
711 include 'COMMON.NAMES'
712 include 'COMMON.CHAIN'
713 include 'COMMON.FFIELD'
714 include 'COMMON.SBRIDGE'
715 include 'COMMON.HEADER'
716 include 'COMMON.CONTROL'
717 include 'COMMON.SAXS'
718 include 'COMMON.DBASE'
719 include 'COMMON.THREAD'
720 include 'COMMON.CONTACTS'
721 include 'COMMON.TORCNSTR'
722 include 'COMMON.TIME1'
723 include 'COMMON.BOUNDS'
725 include 'COMMON.SETUP'
726 include 'COMMON.SHIELD'
727 character*4 sequence(maxres)
729 double precision x(maxvar)
730 character*256 pdbfile
731 character*400 weightcard
732 character*80 weightcard_t,ucase
733 integer itype_pdb(maxres)
734 common /pizda/ itype_pdb
735 logical seq_comp,fail
736 double precision energia(0:n_ene)
737 double precision secprob(3,maxdih_constr)
739 double precision phihel,phibet,sigmahel,sigmabet
740 integer iti,nsi,maxsi
744 integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp
745 double precision sumv
747 C Read PDB structure if applicable
749 if (indpdb.gt.0 .or. pdbref) then
750 read(inp,'(a)') pdbfile
751 if(me.eq.king.or..not.out1file)
752 & write (iout,'(2a)') 'PDB data will be read from file ',
753 & pdbfile(:ilen(pdbfile))
754 open(ipdbin,file=pdbfile,status='old',err=33)
756 33 write (iout,'(a)') 'Error opening PDB file.'
767 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
768 & (crefjlee(j,i+nres),j=1,3)
771 if(me.eq.king.or..not.out1file)
772 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
773 & ' nstart_sup=',nstart_sup
775 itype_pdb(i)=itype(i)
779 nct=nstart_sup+nsup-1
780 call contact(.false.,ncont_ref,icont_ref,co)
783 if(me.eq.king.or..not.out1file)
784 & write(iout,*)'Adding sidechains'
788 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
791 do while (fail.and.nsi.le.maxsi)
792 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
795 if(fail) write(iout,*)'Adding sidechain failed for res ',
796 & i,' after ',nsi,' trials'
801 if (indpdb.eq.0) then
802 C Read sequence if not taken from the pdb file.
804 c print *,'nres=',nres
805 if (iscode.gt.0) then
806 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
808 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
810 C Convert sequence to numeric code
812 itype(i)=rescode(i,sequence(i),iscode)
816 c print '(20i4)',(itype(i),i=1,nres)
819 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
821 if (itype(i).eq.ntyp1) then
825 else if (iabs(itype(i+1)).ne.20) then
827 else if (iabs(itype(i)).ne.20) then
834 if(me.eq.king.or..not.out1file)then
835 write (iout,*) "ITEL"
837 write (iout,*) i,itype(i),itel(i)
839 c print *,'Call Read_Bridge.'
843 cd print *,'NNT=',NNT,' NCT=',NCT
844 call seq2chains(nres,itype,nchain,chain_length,chain_border,
847 chain_border1(2,1)=chain_border(2,1)+1
849 chain_border1(1,i)=chain_border(1,i)-1
850 chain_border1(2,i)=chain_border(2,i)+1
852 if (nchain.gt.1) chain_border1(1,nchain)=chain_border(1,nchain)-1
853 chain_border1(2,nchain)=nres
854 write(iout,*) "nres",nres," nchain",nchain
856 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
857 & chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
859 call chain_symmetry(nchain,nres,itype,chain_border,
860 & chain_length,npermchain,tabpermchain)
862 c write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
865 write(iout,*) "residue permutations"
867 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
870 if (itype(1).eq.ntyp1) nnt=2
871 if (itype(nres).eq.ntyp1) nct=nct-1
872 write (iout,*) "nnt",nnt," nct",nct
875 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
876 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
878 c print*, 'init_dfa_vars finished!'
880 c print*, 'read_dfa_info finished!'
884 if(me.eq.king.or..not.out1file)
885 & write (iout,'(a,i3)') 'nsup=',nsup
887 if (nsup.le.(nct-nnt+1)) then
888 do i=0,nct-nnt+1-nsup
889 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
895 & 'Error - sequences to be superposed do not match.'
898 do i=0,nsup-(nct-nnt+1)
899 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
901 nstart_sup=nstart_sup+i
907 & 'Error - sequences to be superposed do not match.'
910 if (nsup.eq.0) nsup=nct-nnt
911 if (nstart_sup.eq.0) nstart_sup=nnt
912 if (nstart_seq.eq.0) nstart_seq=nnt
913 if(me.eq.king.or..not.out1file)
914 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
915 & ' nstart_seq=',nstart_seq
918 C 8/13/98 Set limits to generating the dihedral angles
923 read (inp,*) ndih_constr
924 write (iout,*) "ndih_constr",ndih_constr
925 if (ndih_constr.gt.0) then
928 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
931 if(me.eq.king.or..not.out1file)then
934 & 'There are',ndih_constr,' restraints on gamma angles.'
936 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
943 phi0(i)=deg2rad*phi0(i)
944 drange(i)=deg2rad*drange(i)
946 C if(me.eq.king.or..not.out1file)
947 C & write (iout,*) 'FTORS',ftors
950 phibound(1,ii) = phi0(i)-drange(i)
951 phibound(2,ii) = phi0(i)+drange(i)
954 if (me.eq.king .or. .not.out1file) then
956 write (iout,'(a)') 'Boundaries in gamma angle sampling:'
958 write (iout,'(a3,i5,2f10.1)')
959 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
964 else if (ndih_constr.lt.0) then
966 call card_concat(weightcard)
967 call reada(weightcard,"PHIHEL",phihel,50.0D0)
968 call reada(weightcard,"PHIBET",phibet,180.0D0)
969 call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
970 call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
971 call reada(weightcard,"WDIHC",wdihc,0.591D0)
972 write (iout,*) "Weight of dihedral angle restraints",wdihc
973 read(inp,'(9x,3f7.3)')
974 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
975 write (iout,*) "The secprob array"
977 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
981 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
982 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
983 ndih_constr=ndih_constr+1
984 idih_constr(ndih_constr)=i
987 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
988 sumv=sumv+vpsipred(j,ndih_constr)
991 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
993 phibound(1,ndih_constr)=phihel*deg2rad
994 phibound(2,ndih_constr)=phibet*deg2rad
995 sdihed(1,ndih_constr)=sigmahel*deg2rad
996 sdihed(2,ndih_constr)=sigmabet*deg2rad
1000 if(me.eq.king.or..not.out1file)then
1003 & 'There are',ndih_constr,
1004 & ' bimodal restraints on gamma angles.'
1006 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
1007 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
1008 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
1009 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
1010 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
1011 & (vpsipred(j,i),j=1,3)
1018 C first setting the theta boundaries to 0 to pi
1019 C this mean that there is no energy penalty for any angle occuring this can be applied
1020 C for generate random conformation but is not implemented in this way
1023 C thetabound(2,i)=pi
1025 C begin reading theta constrains this is quartic constrains allowing to
1026 C have smooth second derivative
1027 if (with_theta_constr) then
1028 C with_theta_constr is keyword allowing for occurance of theta constrains
1029 read (inp,*) ntheta_constr
1030 C ntheta_constr is the number of theta constrains
1031 if (ntheta_constr.gt.0) then
1032 C read (inp,*) ftors
1033 read (inp,*) (itheta_constr(i),theta_constr0(i),
1034 & theta_drange(i),for_thet_constr(i),
1035 & i=1,ntheta_constr)
1036 C the above code reads from 1 to ntheta_constr
1037 C itheta_constr(i) residue i for which is theta_constr
1038 C theta_constr0 the global minimum value
1039 C theta_drange is range for which there is no energy penalty
1040 C for_thet_constr is the force constant for quartic energy penalty
1042 if(me.eq.king.or..not.out1file)then
1044 & 'There are',ntheta_constr,' constraints on phi angles.'
1045 do i=1,ntheta_constr
1046 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1048 & for_thet_constr(i)
1051 do i=1,ntheta_constr
1052 theta_constr0(i)=deg2rad*theta_constr0(i)
1053 theta_drange(i)=deg2rad*theta_drange(i)
1055 C if(me.eq.king.or..not.out1file)
1056 C & write (iout,*) 'FTORS',ftors
1057 C do i=1,ntheta_constr
1058 C ii = itheta_constr(i)
1059 C thetabound(1,ii) = phi0(i)-drange(i)
1060 C thetabound(2,ii) = phi0(i)+drange(i)
1062 endif ! ntheta_constr.gt.0
1063 endif! with_theta_constr
1065 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1066 C write (iout,*) "with_dihed_constr ",with_dihed_constr
1067 c--- Zscore rms -------
1068 if (nz_start.eq.0) nz_start=nnt
1069 if (nz_end.eq.0 .and. nsup.gt.0) then
1071 else if (nz_end.eq.0) then
1074 if(me.eq.king.or..not.out1file)then
1075 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1076 write (iout,*) 'IZ_SC=',iz_sc
1078 c----------------------
1081 if (.not.pdbref) then
1082 call read_angles(inp,*38)
1085 38 write (iout,'(a)') 'Error reading reference structure.'
1087 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1088 stop 'Error reading reference structure'
1090 39 call chainbuild_extconf
1092 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1101 call contact(.true.,ncont_ref,icont_ref,co)
1105 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1107 if (constr_dist.gt.0) call read_dist_constr
1108 c write (iout,*) "After read_dist_constr nhpb",nhpb
1109 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1111 call NMRpeak_partition
1112 if(me.eq.king.or..not.out1file)write (iout,*) 'Contact order:',co
1114 if(me.eq.king.or..not.out1file)
1115 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1118 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1120 if(me.eq.king.or..not.out1file)
1121 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1122 & icont_ref(1,i),' ',
1123 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1126 write (iout,*) "calling read_saxs_consrtr",nsaxs
1127 if (nsaxs.gt.0) call read_saxs_constr
1129 if (constr_homology.gt.0) then
1130 call read_constr_homology
1131 if (indpdb.gt.0 .or. pdbref) then
1134 c(j,i)=crefjlee(j,i)
1135 cref(j,i)=crefjlee(j,i)
1140 write (iout,*) "sc_loc_geom: Array C"
1142 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1143 & (c(j,i+nres),j=1,3)
1145 write (iout,*) "Array Cref"
1147 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1148 & (cref(j,i+nres),j=1,3)
1151 call int_from_cart1(.false.)
1152 call sc_loc_geom(.false.)
1154 thetaref(i)=theta(i)
1159 dc(j,i)=c(j,i+1)-c(j,i)
1160 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1165 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1166 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1171 if (start_from_model) then
1174 read(inp,'(a)',end=332,err=332) pdbfile
1175 if (me.eq.king .or. .not. out1file)
1176 & write (iout,'(a,5x,a)') 'Opening PDB file',
1177 & pdbfile(:ilen(pdbfile))
1178 open(ipdbin,file=pdbfile,status='old',err=336)
1180 336 write (iout,'(a,5x,a)') 'Error opening PDB file',
1181 & pdbfile(:ilen(pdbfile))
1189 if (nres.ge.nres_temp) then
1190 nmodel_start=nmodel_start+1
1191 pdbfiles_chomo(nmodel_start)=pdbfile
1194 chomo(j,i,nmodel_start)=c(j,i)
1198 if (me.eq.king .or. .not. out1file)
1199 & write (iout,'(a,2i5,1x,a)')
1200 & "Different number of residues",nres_temp,nres,
1206 if (nmodel_start.eq.0) then
1207 if (me.eq.king .or. .not. out1file)
1208 & write (iout,'(a)')
1209 & "No valid starting model found START_FROM_MODELS is OFF"
1210 start_from_model=.false.
1212 write (iout,*) "nmodel_start",nmodel_start
1218 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1219 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1220 & modecalc.ne.10) then
1221 C If input structure hasn't been supplied from the PDB file read or generate
1223 if (iranconf.eq.0 .and. .not. extconf .and. .not.
1224 & start_from_model) then
1225 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1226 & write (iout,'(a)') 'Initial geometry will be read in.'
1228 read(inp,'(8f10.5)',end=36,err=36)
1229 & ((c(l,k),l=1,3),k=1,nres),
1230 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1231 if (nnt.gt.1) c(:,nres+1)=c(:,1)
1232 if (nct.lt.nres) c(:,2*nres)=c(:,nres)
1233 c write (iout,*) "Exit READ_CART"
1234 c write (iout,'(8f10.5)')
1235 c & ((c(l,k),l=1,3),k=1,nres),
1236 c & ((c(l,k+nres),l=1,3),k=nnt,nct)
1243 dc(j,i)=c(j,i+1)-c(j,i)
1244 c dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1248 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1250 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1251 c dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1256 c dc_norm(j,i+nres)=0.0d0
1260 call int_from_cart1(.true.)
1261 write (iout,*) "Finish INT_TO_CART"
1262 c write (iout,*) "DC"
1264 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1265 c & (dc(j,i+nres),j=1,3)
1271 call read_angles(inp,*36)
1273 call chainbuild_extconf
1276 36 write (iout,'(a)') 'Error reading angle file.'
1278 call mpi_finalize( MPI_COMM_WORLD,IERR )
1280 stop 'Error reading angle file.'
1282 else if (extconf) then
1283 if (me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1284 & write (iout,'(a)') 'Extended chain initial geometry.'
1286 theta(i)=90d0*deg2rad
1289 phi(i)=180d0*deg2rad
1292 alph(i)=110d0*deg2rad
1295 omeg(i)=-120d0*deg2rad
1296 if (itype(i).le.0) omeg(i)=-omeg(i)
1299 call chainbuild_extconf
1300 else if (.not. start_from_model) then
1301 if(me.eq.king.or..not.out1file)
1302 & write (iout,'(a)') 'Random-generated initial geometry.'
1305 if (me.eq.king .or. fg_rank.eq.0 .and. (
1306 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1310 call gen_rand_conf(itmp,*30)
1312 30 write (iout,*) 'Failed to generate random conformation',
1313 & ', itrial=',itrial
1314 write (*,*) 'Processor:',me,
1315 & ' Failed to generate random conformation',
1325 write (iout,'(a,i3,a)') 'Processor:',me,
1326 & ' error in generating random conformation.'
1327 write (*,'(a,i3,a)') 'Processor:',me,
1328 & ' error in generating random conformation.'
1331 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1336 & ' error in generating random conformation.'
1341 elseif (modecalc.eq.4) then
1342 read (inp,'(a)') intinname
1343 open (intin,file=intinname,status='old',err=333)
1344 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1345 & write (iout,'(a)') 'intinname',intinname
1346 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1348 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1350 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1352 stop 'Error opening angle file.'
1356 C Generate distance constraints, if the PDB structure is to be regularized.
1357 if (nthread.gt.0) then
1358 call read_threadbase
1361 if (me.eq.king .or. .not. out1file)
1363 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1364 write (iout,'(/a,i3,a)')
1365 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1366 write (iout,'(20i4)') (iss(i),i=1,ns)
1368 write(iout,*)"Running with dynamic disulfide-bond formation"
1370 write (iout,'(/a/)') 'Pre-formed links are:'
1376 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1377 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1383 if (ns.gt.0.and.dyn_ss) then
1387 forcon(i-nss)=forcon(i)
1394 dyn_ss_mask(iss(i))=.true.
1399 c write (iout,*) "DC"
1401 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1402 c & (dc(j,i+nres),j=1,3)
1404 if (i2ndstr.gt.0) call secstrp2dihc
1405 c call geom_to_var(nvar,x)
1406 c call etotal(energia(0))
1407 c call enerprint(energia(0))
1408 c call briefout(0,etot)
1410 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1411 cd write (iout,'(a)') 'Variable list:'
1412 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1414 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1415 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1416 & 'Processor',myrank,': end reading molecular data.'
1421 c--------------------------------------------------------------------------
1422 logical function seq_comp(itypea,itypeb,length)
1424 integer length,itypea(length),itypeb(length)
1427 if (itypea(i).ne.itypeb(i)) then
1435 c-----------------------------------------------------------------------------
1436 subroutine read_bridge
1437 C Read information about disulfide bridges.
1439 include 'DIMENSIONS'
1444 include 'COMMON.IOUNITS'
1445 include 'COMMON.GEO'
1446 include 'COMMON.VAR'
1447 include 'COMMON.INTERACT'
1448 include 'COMMON.LOCAL'
1449 include 'COMMON.NAMES'
1450 include 'COMMON.CHAIN'
1451 include 'COMMON.FFIELD'
1452 include 'COMMON.SBRIDGE'
1453 include 'COMMON.HEADER'
1454 include 'COMMON.CONTROL'
1455 include 'COMMON.DBASE'
1456 include 'COMMON.THREAD'
1457 include 'COMMON.TIME1'
1458 include 'COMMON.SETUP'
1460 C Read bridging residues.
1461 read (inp,*) ns,(iss(i),i=1,ns)
1462 c 5/24/2020 Adam: Added a table to translate residue numbers to cysteine
1469 if(me.eq.king.or..not.out1file)
1470 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1471 C Check whether the specified bridging residues are cystines.
1473 if (itype(iss(i)).ne.1) then
1474 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1475 & 'Do you REALLY think that the residue ',
1476 & restyp(itype(iss(i))),i,
1477 & ' can form a disulfide bridge?!!!'
1478 write (*,'(2a,i3,a)')
1479 & 'Do you REALLY think that the residue ',
1480 & restyp(itype(iss(i))),i,
1481 & ' can form a disulfide bridge?!!!'
1483 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1488 C Read preformed bridges.
1490 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1492 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1495 C Check if the residues involved in bridges are in the specified list of
1496 C bridging residues.
1499 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1500 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1501 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1502 & ' contains residues present in other pairs.'
1503 write (*,'(a,i3,a)') 'Disulfide pair',i,
1504 & ' contains residues present in other pairs.'
1506 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1512 if (ihpb(i).eq.iss(j)) goto 10
1514 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1517 if (jhpb(i).eq.iss(j)) goto 20
1519 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1525 ihpb(i)=ihpb(i)+nres
1526 jhpb(i)=jhpb(i)+nres
1532 c----------------------------------------------------------------------------
1533 subroutine read_x(kanal,*)
1535 include 'DIMENSIONS'
1536 include 'COMMON.GEO'
1537 include 'COMMON.VAR'
1538 include 'COMMON.CHAIN'
1539 include 'COMMON.IOUNITS'
1540 include 'COMMON.CONTROL'
1541 include 'COMMON.LOCAL'
1542 include 'COMMON.INTERACT'
1543 integer i,j,k,l,kanal
1544 c Read coordinates from input
1546 read(kanal,'(8f10.5)',end=10,err=10)
1547 & ((c(l,k),l=1,3),k=1,nres),
1548 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1551 c(j,2*nres)=c(j,nres)
1553 call int_from_cart1(.false.)
1556 dc(j,i)=c(j,i+1)-c(j,i)
1557 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1561 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1563 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1564 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1572 c----------------------------------------------------------------------------
1573 subroutine read_threadbase
1575 include 'DIMENSIONS'
1576 include 'COMMON.IOUNITS'
1577 include 'COMMON.GEO'
1578 include 'COMMON.VAR'
1579 include 'COMMON.INTERACT'
1580 include 'COMMON.LOCAL'
1581 include 'COMMON.NAMES'
1582 include 'COMMON.CHAIN'
1583 include 'COMMON.FFIELD'
1584 include 'COMMON.SBRIDGE'
1585 include 'COMMON.HEADER'
1586 include 'COMMON.CONTROL'
1587 include 'COMMON.DBASE'
1588 include 'COMMON.THREAD'
1589 include 'COMMON.TIME1'
1591 double precision dist
1592 C Read pattern database for threading.
1593 read (icbase,*) nseq
1595 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1596 & nres_base(2,i),nres_base(3,i)
1597 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1599 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1600 c & nres_base(2,i),nres_base(3,i)
1601 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1605 if (weidis.eq.0.0D0) weidis=0.1D0
1614 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1615 write (iout,'(a,i5)') 'nexcl: ',nexcl
1616 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1619 c------------------------------------------------------------------------------
1620 subroutine setup_var
1622 include 'DIMENSIONS'
1623 include 'COMMON.IOUNITS'
1624 include 'COMMON.GEO'
1625 include 'COMMON.VAR'
1626 include 'COMMON.INTERACT'
1627 include 'COMMON.LOCAL'
1628 include 'COMMON.NAMES'
1629 include 'COMMON.CHAIN'
1630 include 'COMMON.FFIELD'
1631 include 'COMMON.SBRIDGE'
1632 include 'COMMON.HEADER'
1633 include 'COMMON.CONTROL'
1634 include 'COMMON.DBASE'
1635 include 'COMMON.THREAD'
1636 include 'COMMON.TIME1'
1638 C Set up variable list.
1644 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1646 ialph(i,1)=nvar+nside
1650 if (indphi.gt.0) then
1652 else if (indback.gt.0) then
1659 c----------------------------------------------------------------------------
1660 subroutine gen_dist_constr
1661 C Generate CA distance constraints.
1663 include 'DIMENSIONS'
1664 include 'COMMON.IOUNITS'
1665 include 'COMMON.GEO'
1666 include 'COMMON.VAR'
1667 include 'COMMON.INTERACT'
1668 include 'COMMON.LOCAL'
1669 include 'COMMON.NAMES'
1670 include 'COMMON.CHAIN'
1671 include 'COMMON.FFIELD'
1672 include 'COMMON.SBRIDGE'
1673 include 'COMMON.HEADER'
1674 include 'COMMON.CONTROL'
1675 include 'COMMON.DBASE'
1676 include 'COMMON.THREAD'
1677 include 'COMMON.SPLITELE'
1678 include 'COMMON.TIME1'
1679 integer i,j,itype_pdb(maxres)
1680 common /pizda/ itype_pdb
1682 double precision dist
1684 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1685 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1686 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1688 do i=nstart_sup,nstart_sup+nsup-1
1689 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1690 cd & ' seq_pdb', restyp(itype_pdb(i))
1691 do j=i+2,nstart_sup+nsup-1
1692 c 5/24/2020 Adam: Cutoff included to reduce array size
1694 if (dd.gt.r_cut_int) cycle
1696 ihpb(nhpb)=i+nstart_seq-nstart_sup
1697 jhpb(nhpb)=j+nstart_seq-nstart_sup
1702 cd write (iout,'(a)') 'Distance constraints:'
1707 cd if (ii.gt.nres) then
1712 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1713 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1714 cd & dhpb(i),forcon(i)
1718 c----------------------------------------------------------------------------
1721 include 'DIMENSIONS'
1722 include 'COMMON.MAP'
1723 include 'COMMON.IOUNITS'
1725 character*3 angid(4) /'THE','PHI','ALP','OME'/
1726 character*80 mapcard,ucase
1728 read (inp,'(a)') mapcard
1729 mapcard=ucase(mapcard)
1730 if (index(mapcard,'PHI').gt.0) then
1732 else if (index(mapcard,'THE').gt.0) then
1734 else if (index(mapcard,'ALP').gt.0) then
1736 else if (index(mapcard,'OME').gt.0) then
1739 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1740 stop 'Error - illegal variable spec in MAP card.'
1742 call readi (mapcard,'RES1',res1(imap),0)
1743 call readi (mapcard,'RES2',res2(imap),0)
1744 if (res1(imap).eq.0) then
1745 res1(imap)=res2(imap)
1746 else if (res2(imap).eq.0) then
1747 res2(imap)=res1(imap)
1749 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1751 & 'Error - illegal definition of variable group in MAP.'
1752 stop 'Error - illegal definition of variable group in MAP.'
1754 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1755 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1756 call readi(mapcard,'NSTEP',nstep(imap),0)
1757 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1759 & 'Illegal boundary and/or step size specification in MAP.'
1760 stop 'Illegal boundary and/or step size specification in MAP.'
1765 c----------------------------------------------------------------------------
1768 include 'DIMENSIONS'
1769 include 'COMMON.IOUNITS'
1770 include 'COMMON.GEO'
1771 include 'COMMON.CSA'
1772 include 'COMMON.BANK'
1773 include 'COMMON.CONTROL'
1775 character*620 mcmcard
1776 call card_concat(mcmcard)
1778 call readi(mcmcard,'NCONF',nconf,50)
1779 call readi(mcmcard,'NADD',nadd,0)
1780 call readi(mcmcard,'JSTART',jstart,1)
1781 call readi(mcmcard,'JEND',jend,1)
1782 call readi(mcmcard,'NSTMAX',nstmax,500000)
1783 call readi(mcmcard,'N0',n0,1)
1784 call readi(mcmcard,'N1',n1,6)
1785 call readi(mcmcard,'N2',n2,4)
1786 call readi(mcmcard,'N3',n3,0)
1787 call readi(mcmcard,'N4',n4,0)
1788 call readi(mcmcard,'N5',n5,0)
1789 call readi(mcmcard,'N6',n6,10)
1790 call readi(mcmcard,'N7',n7,0)
1791 call readi(mcmcard,'N8',n8,0)
1792 call readi(mcmcard,'N9',n9,0)
1793 call readi(mcmcard,'N14',n14,0)
1794 call readi(mcmcard,'N15',n15,0)
1795 call readi(mcmcard,'N16',n16,0)
1796 call readi(mcmcard,'N17',n17,0)
1797 call readi(mcmcard,'N18',n18,0)
1799 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1801 call readi(mcmcard,'NDIFF',ndiff,2)
1802 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1803 call readi(mcmcard,'IS1',is1,1)
1804 call readi(mcmcard,'IS2',is2,8)
1805 call readi(mcmcard,'NRAN0',nran0,4)
1806 call readi(mcmcard,'NRAN1',nran1,2)
1807 call readi(mcmcard,'IRR',irr,1)
1808 call readi(mcmcard,'NSEED',nseed,20)
1809 call readi(mcmcard,'NTOTAL',ntotal,10000)
1810 call reada(mcmcard,'CUT1',cut1,2.0d0)
1811 call reada(mcmcard,'CUT2',cut2,5.0d0)
1812 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1813 call readi(mcmcard,'ICMAX',icmax,3)
1814 call readi(mcmcard,'IRESTART',irestart,0)
1815 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1818 call reada(mcmcard,'DELE',dele,20.0d0)
1819 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1820 call readi(mcmcard,'IREF',iref,0)
1821 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1822 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1823 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1824 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1825 write (iout,*) "NCONF_IN",nconf_in
1828 c----------------------------------------------------------------------------
1831 include 'DIMENSIONS'
1832 include 'COMMON.MCM'
1833 include 'COMMON.MCE'
1834 include 'COMMON.IOUNITS'
1836 character*320 mcmcard
1838 call card_concat(mcmcard)
1839 call readi(mcmcard,'MAXACC',maxacc,100)
1840 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1841 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1842 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1843 call readi(mcmcard,'MAXREPM',maxrepm,200)
1844 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1845 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1846 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1847 call reada(mcmcard,'E_UP',e_up,5.0D0)
1848 call reada(mcmcard,'DELTE',delte,0.1D0)
1849 call readi(mcmcard,'NSWEEP',nsweep,5)
1850 call readi(mcmcard,'NSTEPH',nsteph,0)
1851 call readi(mcmcard,'NSTEPC',nstepc,0)
1852 call reada(mcmcard,'TMIN',tmin,298.0D0)
1853 call reada(mcmcard,'TMAX',tmax,298.0D0)
1854 call readi(mcmcard,'NWINDOW',nwindow,0)
1855 call readi(mcmcard,'PRINT_MC',print_mc,0)
1856 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1857 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1858 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1859 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1860 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1861 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1862 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1863 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1864 if (nwindow.gt.0) then
1865 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1867 winlen(i)=winend(i)-winstart(i)+1
1870 if (tmax.lt.tmin) tmax=tmin
1871 if (tmax.eq.tmin) then
1875 if (nstepc.gt.0 .and. nsteph.gt.0) then
1876 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1877 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1879 C Probabilities of different move types
1880 sumpro_type(0)=0.0D0
1881 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1882 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1883 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1884 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1885 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1886 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1887 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1889 print *,'i',i,' sumprotype',sumpro_type(i)
1890 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1891 print *,'i',i,' sumprotype',sumpro_type(i)
1895 c----------------------------------------------------------------------------
1896 subroutine read_minim
1898 include 'DIMENSIONS'
1899 include 'COMMON.MINIM'
1900 include 'COMMON.IOUNITS'
1902 character*320 minimcard
1903 call card_concat(minimcard)
1904 call readi(minimcard,'MAXMIN',maxmin,2000)
1905 call readi(minimcard,'MAXFUN',maxfun,5000)
1906 call readi(minimcard,'MINMIN',minmin,maxmin)
1907 call readi(minimcard,'MINFUN',minfun,maxmin)
1908 call reada(minimcard,'TOLF',tolf,1.0D-2)
1909 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1910 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1911 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1912 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1913 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1914 & 'Options in energy minimization:'
1915 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1916 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1917 & 'MinMin:',MinMin,' MinFun:',MinFun,
1918 & ' TolF:',TolF,' RTolF:',RTolF
1921 c----------------------------------------------------------------------------
1922 subroutine read_angles(kanal,*)
1924 include 'DIMENSIONS'
1925 include 'COMMON.GEO'
1926 include 'COMMON.VAR'
1927 include 'COMMON.CHAIN'
1928 include 'COMMON.IOUNITS'
1929 include 'COMMON.CONTROL'
1931 c Read angles from input
1933 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1934 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1935 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1936 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1939 c 9/7/01 avoid 180 deg valence angle
1940 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1942 theta(i)=deg2rad*theta(i)
1943 phi(i)=deg2rad*phi(i)
1944 alph(i)=deg2rad*alph(i)
1945 omeg(i)=deg2rad*omeg(i)
1950 c----------------------------------------------------------------------------
1951 subroutine reada(rekord,lancuch,wartosc,default)
1953 character*(*) rekord,lancuch
1954 double precision wartosc,default
1957 iread=index(rekord,lancuch)
1958 if (iread.eq.0) then
1962 iread=iread+ilen(lancuch)+1
1963 read (rekord(iread:),*,err=10,end=10) wartosc
1968 c----------------------------------------------------------------------------
1969 subroutine readi(rekord,lancuch,wartosc,default)
1971 character*(*) rekord,lancuch
1972 integer wartosc,default
1975 iread=index(rekord,lancuch)
1976 if (iread.eq.0) then
1980 iread=iread+ilen(lancuch)+1
1981 read (rekord(iread:),*,err=10,end=10) wartosc
1986 c----------------------------------------------------------------------------
1987 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1990 integer tablica(dim),default
1991 character*(*) rekord,lancuch
1998 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1999 if (iread.eq.0) return
2000 iread=iread+ilen(lancuch)+1
2001 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2004 c----------------------------------------------------------------------------
2005 subroutine multreada(rekord,lancuch,tablica,dim,default)
2008 double precision tablica(dim),default
2009 character*(*) rekord,lancuch
2016 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2017 if (iread.eq.0) return
2018 iread=iread+ilen(lancuch)+1
2019 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2022 c----------------------------------------------------------------------------
2023 subroutine openunits
2025 include 'DIMENSIONS'
2029 character*16 form,nodename
2032 include 'COMMON.SETUP'
2033 include 'COMMON.IOUNITS'
2035 include 'COMMON.CONTROL'
2036 integer lenpre,lenpot,ilen,lentmp,npos
2038 character*3 out1file_text,ucase
2043 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2044 call getenv_loc("PREFIX",prefix)
2046 call getenv_loc("POT",pot)
2047 call getenv_loc("DIRTMP",tmpdir)
2048 call getenv_loc("CURDIR",curdir)
2049 call getenv_loc("OUT1FILE",out1file_text)
2050 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2051 out1file_text=ucase(out1file_text)
2052 if (out1file_text(1:1).eq."Y") then
2055 out1file=fg_rank.gt.0
2060 if (lentmp.gt.0) then
2061 write (*,'(80(1h!))')
2062 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2063 write (*,'(80(1h!))')
2064 write (*,*)"All output files will be on node /tmp directory."
2066 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2067 if (me.eq.king) then
2068 write (*,*) "The master node is ",nodename
2069 else if (fg_rank.eq.0) then
2070 write (*,*) "I am the CG slave node ",nodename
2072 write (*,*) "I am the FG slave node ",nodename
2075 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2076 lenpre = lentmp+lenpre+1
2078 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2079 C Get the names and open the input files
2080 #if defined(WINIFL) || defined(WINPGI)
2081 open(1,file=pref_orig(:ilen(pref_orig))//
2082 & '.inp',status='old',readonly,shared)
2083 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2084 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2085 C Get parameter filenames and open the parameter files.
2086 call getenv_loc('BONDPAR',bondname)
2087 open (ibond,file=bondname,status='old',readonly,shared)
2088 call getenv_loc('THETPAR',thetname)
2089 open (ithep,file=thetname,status='old',readonly,shared)
2090 call getenv_loc('ROTPAR',rotname)
2091 open (irotam,file=rotname,status='old',readonly,shared)
2092 call getenv_loc('TORPAR',torname)
2093 open (itorp,file=torname,status='old',readonly,shared)
2094 call getenv_loc('TORDPAR',tordname)
2095 open (itordp,file=tordname,status='old',readonly,shared)
2096 call getenv_loc('FOURIER',fouriername)
2097 open (ifourier,file=fouriername,status='old',readonly,shared)
2098 call getenv_loc('ELEPAR',elename)
2099 open (ielep,file=elename,status='old',readonly,shared)
2100 call getenv_loc('SIDEPAR',sidename)
2101 open (isidep,file=sidename,status='old',readonly,shared)
2102 call getenv_loc('LIPTRANPAR',liptranname)
2103 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2104 #elif (defined CRAY) || (defined AIX)
2105 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2107 c print *,"Processor",myrank," opened file 1"
2108 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2109 c print *,"Processor",myrank," opened file 9"
2110 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2111 C Get parameter filenames and open the parameter files.
2112 call getenv_loc('BONDPAR',bondname)
2113 open (ibond,file=bondname,status='old',action='read')
2114 c print *,"Processor",myrank," opened file IBOND"
2115 call getenv_loc('THETPAR',thetname)
2116 open (ithep,file=thetname,status='old',action='read')
2118 call getenv_loc('THETPARPDB',thetname_pdb)
2119 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2121 c print *,"Processor",myrank," opened file ITHEP"
2122 call getenv_loc('ROTPAR',rotname)
2123 open (irotam,file=rotname,status='old',action='read')
2125 call getenv_loc('ROTPARPDB',rotname_pdb)
2126 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2128 c print *,"Processor",myrank," opened file IROTAM"
2129 call getenv_loc('TORPAR',torname)
2130 open (itorp,file=torname,status='old',action='read')
2131 c print *,"Processor",myrank," opened file ITORP"
2132 call getenv_loc('TORDPAR',tordname)
2133 open (itordp,file=tordname,status='old',action='read')
2134 c print *,"Processor",myrank," opened file ITORDP"
2135 call getenv_loc('SCCORPAR',sccorname)
2136 open (isccor,file=sccorname,status='old',action='read')
2137 c print *,"Processor",myrank," opened file ISCCOR"
2138 call getenv_loc('FOURIER',fouriername)
2139 open (ifourier,file=fouriername,status='old',action='read')
2140 c print *,"Processor",myrank," opened file IFOURIER"
2141 call getenv_loc('ELEPAR',elename)
2142 open (ielep,file=elename,status='old',action='read')
2143 c print *,"Processor",myrank," opened file IELEP"
2144 call getenv_loc('SIDEPAR',sidename)
2145 open (isidep,file=sidename,status='old',action='read')
2146 call getenv_loc('LIPTRANPAR',liptranname)
2147 open (iliptranpar,file=liptranname,status='old',action='read')
2148 c print *,"Processor",myrank," opened file ISIDEP"
2149 c print *,"Processor",myrank," opened parameter files"
2151 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2152 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2153 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2154 C Get parameter filenames and open the parameter files.
2155 call getenv_loc('BONDPAR',bondname)
2156 open (ibond,file=bondname,status='old')
2157 call getenv_loc('THETPAR',thetname)
2158 open (ithep,file=thetname,status='old')
2160 call getenv_loc('THETPARPDB',thetname_pdb)
2161 open (ithep_pdb,file=thetname_pdb,status='old')
2163 call getenv_loc('ROTPAR',rotname)
2164 open (irotam,file=rotname,status='old')
2166 call getenv_loc('ROTPARPDB',rotname_pdb)
2167 open (irotam_pdb,file=rotname_pdb,status='old')
2169 call getenv_loc('TORPAR',torname)
2170 open (itorp,file=torname,status='old')
2171 call getenv_loc('TORDPAR',tordname)
2172 open (itordp,file=tordname,status='old')
2173 call getenv_loc('SCCORPAR',sccorname)
2174 open (isccor,file=sccorname,status='old')
2175 call getenv_loc('FOURIER',fouriername)
2176 open (ifourier,file=fouriername,status='old')
2177 call getenv_loc('ELEPAR',elename)
2178 open (ielep,file=elename,status='old')
2179 call getenv_loc('SIDEPAR',sidename)
2180 open (isidep,file=sidename,status='old')
2181 call getenv_loc('LIPTRANPAR',liptranname)
2182 open (iliptranpar,file=liptranname,status='old')
2184 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2186 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2187 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2188 C Get parameter filenames and open the parameter files.
2189 call getenv_loc('BONDPAR',bondname)
2190 open (ibond,file=bondname,status='old',readonly)
2191 call getenv_loc('THETPAR',thetname)
2192 open (ithep,file=thetname,status='old',readonly)
2193 call getenv_loc('ROTPAR',rotname)
2194 open (irotam,file=rotname,status='old',readonly)
2195 call getenv_loc('TORPAR',torname)
2196 open (itorp,file=torname,status='old',readonly)
2197 call getenv_loc('TORDPAR',tordname)
2198 open (itordp,file=tordname,status='old',readonly)
2199 call getenv_loc('SCCORPAR',sccorname)
2200 open (isccor,file=sccorname,status='old',readonly)
2202 call getenv_loc('THETPARPDB',thetname_pdb)
2203 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2205 call getenv_loc('FOURIER',fouriername)
2206 open (ifourier,file=fouriername,status='old',readonly)
2207 call getenv_loc('ELEPAR',elename)
2208 open (ielep,file=elename,status='old',readonly)
2209 call getenv_loc('SIDEPAR',sidename)
2210 open (isidep,file=sidename,status='old',readonly)
2211 call getenv_loc('LIPTRANPAR',liptranname)
2212 open (iliptranpar,file=liptranname,status='old',action='read')
2214 call getenv_loc('ROTPARPDB',rotname_pdb)
2215 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2219 call getenv_loc('TUBEPAR',tubename)
2220 #if defined(WINIFL) || defined(WINPGI)
2221 open (itube,file=tubename,status='old',readonly,shared)
2222 #elif (defined CRAY) || (defined AIX)
2223 open (itube,file=tubename,status='old',action='read')
2225 open (itube,file=tubename,status='old')
2227 open (itube,file=tubename,status='old',readonly)
2232 C 8/9/01 In the newest version SCp interaction constants are read from a file
2233 C Use -DOLDSCP to use hard-coded constants instead.
2235 call getenv_loc('SCPPAR',scpname)
2236 #if defined(WINIFL) || defined(WINPGI)
2237 open (iscpp,file=scpname,status='old',readonly,shared)
2238 #elif (defined CRAY) || (defined AIX)
2239 open (iscpp,file=scpname,status='old',action='read')
2241 open (iscpp,file=scpname,status='old')
2243 open (iscpp,file=scpname,status='old',readonly)
2246 call getenv_loc('PATTERN',patname)
2247 #if defined(WINIFL) || defined(WINPGI)
2248 open (icbase,file=patname,status='old',readonly,shared)
2249 #elif (defined CRAY) || (defined AIX)
2250 open (icbase,file=patname,status='old',action='read')
2252 open (icbase,file=patname,status='old')
2254 open (icbase,file=patname,status='old',readonly)
2257 C Open output file only for CG processes
2258 c print *,"Processor",myrank," fg_rank",fg_rank
2259 if (fg_rank.eq.0) then
2261 if (nodes.eq.1) then
2264 npos = dlog10(dfloat(nodes-1))+1
2266 if (npos.lt.3) npos=3
2267 write (liczba,'(i1)') npos
2268 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2270 write (liczba,form) me
2271 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2272 & liczba(:ilen(liczba))
2273 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2275 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2277 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2278 & liczba(:ilen(liczba))//'.mol2'
2279 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2280 & liczba(:ilen(liczba))//'.stat'
2282 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2283 & //liczba(:ilen(liczba))//'.stat')
2284 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2287 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2288 & liczba(:ilen(liczba))//'.const'
2293 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2294 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2295 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2296 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2297 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2299 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2301 rest2name=prefix(:ilen(prefix))//'.rst'
2303 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2306 #if defined(AIX) || defined(PGI) || defined(CRAY)
2307 if (me.eq.king .or. .not. out1file) then
2308 open(iout,file=outname,status='unknown')
2310 open(iout,file="/dev/null",status="unknown")
2314 if (fg_rank.gt.0) then
2315 write (liczba,'(i3.3)') myrank/nfgtasks
2316 write (ll,'(bz,i3.3)') fg_rank
2317 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2323 open(igeom,file=intname,status='unknown',position='append')
2324 open(ipdb,file=pdbname,status='unknown')
2325 open(imol2,file=mol2name,status='unknown')
2326 open(istat,file=statname,status='unknown',position='append')
2328 c1out open(iout,file=outname,status='unknown')
2331 if (me.eq.king .or. .not.out1file)
2332 & open(iout,file=outname,status='unknown')
2334 if (fg_rank.gt.0) then
2335 write (liczba,'(i3.3)') myrank/nfgtasks
2336 write (ll,'(bz,i3.3)') fg_rank
2337 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2342 open(igeom,file=intname,status='unknown',access='append')
2343 open(ipdb,file=pdbname,status='unknown')
2344 open(imol2,file=mol2name,status='unknown')
2345 open(istat,file=statname,status='unknown',access='append')
2347 c1out open(iout,file=outname,status='unknown')
2350 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2351 csa_seed=prefix(:lenpre)//'.CSA.seed'
2352 csa_history=prefix(:lenpre)//'.CSA.history'
2353 csa_bank=prefix(:lenpre)//'.CSA.bank'
2354 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2355 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2356 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2357 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2358 csa_int=prefix(:lenpre)//'.int'
2359 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2360 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2361 csa_in=prefix(:lenpre)//'.CSA.in'
2362 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2365 write (iout,'(80(1h-))')
2366 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2367 write (iout,'(80(1h-))')
2368 write (iout,*) "Input file : ",
2369 & pref_orig(:ilen(pref_orig))//'.inp'
2370 write (iout,*) "Output file : ",
2371 & outname(:ilen(outname))
2373 write (iout,*) "Sidechain potential file : ",
2374 & sidename(:ilen(sidename))
2376 write (iout,*) "SCp potential file : ",
2377 & scpname(:ilen(scpname))
2379 write (iout,*) "Electrostatic potential file : ",
2380 & elename(:ilen(elename))
2381 write (iout,*) "Cumulant coefficient file : ",
2382 & fouriername(:ilen(fouriername))
2383 write (iout,*) "Torsional parameter file : ",
2384 & torname(:ilen(torname))
2385 write (iout,*) "Double torsional parameter file : ",
2386 & tordname(:ilen(tordname))
2387 write (iout,*) "SCCOR parameter file : ",
2388 & sccorname(:ilen(sccorname))
2389 write (iout,*) "Bond & inertia constant file : ",
2390 & bondname(:ilen(bondname))
2391 write (iout,*) "Bending-torsion parameter file : ",
2392 & thetname(:ilen(thetname))
2393 write (iout,*) "Rotamer parameter file : ",
2394 & rotname(:ilen(rotname))
2395 write (iout,*) "Thetpdb parameter file : ",
2396 & thetname_pdb(:ilen(thetname_pdb))
2397 write (iout,*) "Threading database : ",
2398 & patname(:ilen(patname))
2400 &write (iout,*)" DIRTMP : ",
2402 write (iout,'(80(1h-))')
2406 c----------------------------------------------------------------------------
2407 subroutine card_concat(card)
2408 implicit real*8 (a-h,o-z)
2409 include 'DIMENSIONS'
2410 include 'COMMON.IOUNITS'
2412 character*80 karta,ucase
2414 read (inp,'(a)') karta
2417 do while (karta(80:80).eq.'&')
2418 card=card(:ilen(card)+1)//karta(:79)
2419 read (inp,'(a)') karta
2422 card=card(:ilen(card)+1)//karta
2425 c------------------------------------------------------------------------------
2428 include 'DIMENSIONS'
2429 include 'COMMON.CHAIN'
2430 include 'COMMON.IOUNITS'
2431 include 'COMMON.CONTROL'
2433 include 'COMMON.QRESTR'
2435 include 'COMMON.LAGRANGE.5diag'
2437 include 'COMMON.LAGRANGE'
2440 open(irest2,file=rest2name,status='unknown')
2441 read(irest2,*) totT,EK,potE,totE,t_bath
2444 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2447 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2450 read (irest2,*) iset
2455 c------------------------------------------------------------------------------
2456 subroutine read_fragments
2458 include 'DIMENSIONS'
2462 include 'COMMON.SETUP'
2463 include 'COMMON.CHAIN'
2464 include 'COMMON.IOUNITS'
2466 include 'COMMON.QRESTR'
2467 include 'COMMON.CONTROL'
2469 read(inp,*) nset,nfrag,npair,nfrag_back
2470 loc_qlike=(nfrag_back.lt.0)
2471 nfrag_back=iabs(nfrag_back)
2472 if(me.eq.king.or..not.out1file)
2473 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2474 & " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2476 read(inp,*) mset(iset)
2478 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2480 if(me.eq.king.or..not.out1file)
2481 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2482 & ifrag(2,i,iset), qinfrag(i,iset)
2485 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2487 if(me.eq.king.or..not.out1file)
2488 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2489 & ipair(2,i,iset), qinpair(i,iset)
2493 read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2494 & wfrag_back(2,i,iset),qin_back(2,i,iset),
2495 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2496 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2497 if(me.eq.king.or..not.out1file)
2498 & write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2499 & wfrag_back(2,i,iset),qin_back(3,i,iset),
2500 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2501 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2505 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2506 & wfrag_back(3,i,iset),
2507 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2508 if(me.eq.king.or..not.out1file)
2509 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2510 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2516 C---------------------------------------------------------------------------
2517 subroutine read_afminp
2519 include 'DIMENSIONS'
2523 include 'COMMON.SETUP'
2524 include 'COMMON.CONTROL'
2525 include 'COMMON.CHAIN'
2526 include 'COMMON.IOUNITS'
2527 include 'COMMON.SBRIDGE'
2528 character*320 afmcard
2530 c print *, "wchodze"
2531 call card_concat(afmcard)
2532 call readi(afmcard,"BEG",afmbeg,0)
2533 call readi(afmcard,"END",afmend,0)
2534 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2535 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2536 c print *,'FORCE=' ,forceAFMconst
2537 CCCC NOW PROPERTIES FOR AFM
2540 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2542 distafminit=dsqrt(distafminit)
2543 c print *,'initdist',distafminit
2546 c-------------------------------------------------------------------------------
2547 subroutine read_saxs_constr
2549 include 'DIMENSIONS'
2553 include 'COMMON.SETUP'
2554 include 'COMMON.CONTROL'
2555 include 'COMMON.SAXS'
2556 include 'COMMON.CHAIN'
2557 include 'COMMON.IOUNITS'
2558 include 'COMMON.SBRIDGE'
2559 double precision cm(3),cnorm
2562 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2564 if (saxs_mode.eq.0) then
2565 c SAXS distance distribution
2567 read(inp,*) distsaxs(i),Psaxs(i)
2571 Cnorm = Cnorm + Psaxs(i)
2573 write (iout,*) "Cnorm",Cnorm
2575 Psaxs(i)=Psaxs(i)/Cnorm
2577 write (iout,*) "Normalized distance distribution from SAXS"
2579 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2583 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2585 write (iout,*) "Wsaxs0",Wsaxs0
2589 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2596 cm(j)=cm(j)+Csaxs(j,i)
2604 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2607 write (iout,*) "SAXS sphere coordinates"
2609 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2615 c-------------------------------------------------------------------------------
2616 subroutine read_dist_constr
2618 include 'DIMENSIONS'
2622 include 'COMMON.SETUP'
2623 include 'COMMON.CONTROL'
2624 include 'COMMON.CHAIN'
2625 include 'COMMON.IOUNITS'
2626 include 'COMMON.SBRIDGE'
2627 include 'COMMON.INTERACT'
2628 integer i,j,k,ii,jj,itemp,link_type,iiend,jjend,kk
2629 integer nfrag_,npair_,ndist_,ifrag_(2,100),ipair_(2,1000)
2630 double precision wfrag_(100),wpair_(1000)
2631 double precision ddjk,dist,dist_cut,fordepthmax
2632 character*5000 controlcard
2633 logical normalize,next
2635 double precision scal_bfac
2636 double precision xlink(4,0:4) /
2638 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
2639 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
2640 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
2641 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
2642 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
2643 c print *, "WCHODZE"
2644 write (iout,*) "Calling read_dist_constr"
2645 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2647 restr_on_coord=.false.
2656 call card_concat(controlcard)
2657 next = index(controlcard,"NEXT").gt.0
2658 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2659 write (iout,*) "restr_type",restr_type
2660 call readi(controlcard,"NFRAG",nfrag_,0)
2661 call readi(controlcard,"NFRAG",nfrag_,0)
2662 call readi(controlcard,"NPAIR",npair_,0)
2663 call readi(controlcard,"NDIST",ndist_,0)
2664 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2665 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2666 if (restr_type.eq.10)
2667 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2668 if (restr_type.eq.12)
2669 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2670 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2671 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2672 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2673 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2674 normalize = index(controlcard,"NORMALIZE").gt.0
2675 write (iout,*) "WBOLTZD",wboltzd
2676 write (iout,*) "SCAL_PEAK",scal_peak
2677 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2678 write (iout,*) "IFRAG"
2680 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2682 write (iout,*) "IPAIR"
2684 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2686 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2688 & "Distance restraints as generated from reference structure"
2690 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2691 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2692 & ifrag_(2,i)=nstart_sup+nsup-1
2693 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2695 if (wfrag_(i).eq.0.0d0) cycle
2696 do j=ifrag_(1,i),ifrag_(2,i)-1
2697 do k=j+1,ifrag_(2,i)
2698 c write (iout,*) "j",j," k",k
2700 if (restr_type.eq.1) then
2706 forcon(nhpb)=wfrag_(i)
2707 else if (constr_dist.eq.2) then
2708 if (ddjk.le.dist_cut) then
2714 forcon(nhpb)=wfrag_(i)
2716 else if (restr_type.eq.3) then
2722 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2725 if (.not.out1file .or. me.eq.king)
2726 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2727 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2729 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2730 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2736 if (wpair_(i).eq.0.0d0) cycle
2744 do j=ifrag_(1,ii),ifrag_(2,ii)
2745 do k=ifrag_(1,jj),ifrag_(2,jj)
2747 if (restr_type.eq.1) then
2753 forcon(nhpb)=wpair_(i)
2754 else if (constr_dist.eq.2) then
2755 if (ddjk.le.dist_cut) then
2761 forcon(nhpb)=wpair_(i)
2763 else if (restr_type.eq.3) then
2769 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2772 if (.not.out1file .or. me.eq.king)
2773 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2774 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2776 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2777 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2784 write (iout,*) "Distance restraints as read from input"
2786 if (restr_type.eq.12) then
2787 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2788 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2789 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2790 & fordepth_peak(nhpb_peak+1),npeak
2791 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2792 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2793 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2794 c & fordepth_peak(nhpb_peak+1),npeak
2795 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2796 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2797 nhpb_peak=nhpb_peak+1
2798 irestr_type_peak(nhpb_peak)=12
2799 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2802 if (.not.out1file .or. me.eq.king)
2803 & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2804 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2805 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2806 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2807 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2809 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2810 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2811 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2812 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2813 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2815 if (ibecarb_peak(nhpb_peak).eq.3) then
2816 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2817 else if (ibecarb_peak(nhpb_peak).eq.2) then
2818 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2819 else if (ibecarb_peak(nhpb_peak).eq.1) then
2820 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2821 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2823 else if (restr_type.eq.11) then
2824 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2825 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2826 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2827 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2829 irestr_type(nhpb)=11
2831 if (.not.out1file .or. me.eq.king)
2832 & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2833 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2834 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2836 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2837 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2838 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2840 c if (ibecarb(nhpb).gt.0) then
2841 c ihpb(nhpb)=ihpb(nhpb)+nres
2842 c jhpb(nhpb)=jhpb(nhpb)+nres
2844 if (ibecarb(nhpb).eq.3) then
2845 ihpb(nhpb)=ihpb(nhpb)+nres
2846 else if (ibecarb(nhpb).eq.2) then
2847 ihpb(nhpb)=ihpb(nhpb)+nres
2848 else if (ibecarb(nhpb).eq.1) then
2849 ihpb(nhpb)=ihpb(nhpb)+nres
2850 jhpb(nhpb)=jhpb(nhpb)+nres
2852 else if (restr_type.eq.10) then
2853 c Cross-lonk Markov-like potential
2854 call card_concat(controlcard)
2855 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2856 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2858 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2859 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2860 if (index(controlcard,"ZL").gt.0) then
2862 else if (index(controlcard,"ADH").gt.0) then
2864 else if (index(controlcard,"PDH").gt.0) then
2866 else if (index(controlcard,"DSS").gt.0) then
2871 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2872 & xlink(1,link_type))
2873 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2874 & xlink(2,link_type))
2875 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2876 & xlink(3,link_type))
2877 call reada(controlcard,"SIGMA",forcon(nhpb+1),
2878 & xlink(4,link_type))
2879 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2880 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2881 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2882 if (forcon(nhpb+1).le.0.0d0 .or.
2883 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2885 irestr_type(nhpb)=10
2886 if (ibecarb(nhpb).eq.3) then
2887 jhpb(nhpb)=jhpb(nhpb)+nres
2888 else if (ibecarb(nhpb).eq.2) then
2889 ihpb(nhpb)=ihpb(nhpb)+nres
2890 else if (ibecarb(nhpb).eq.1) then
2891 ihpb(nhpb)=ihpb(nhpb)+nres
2892 jhpb(nhpb)=jhpb(nhpb)+nres
2895 if (.not.out1file .or. me.eq.king)
2896 & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2897 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2898 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2901 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2902 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2903 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2908 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2909 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2910 if (forcon(nhpb+1).gt.0.0d0) then
2912 if (dhpb1(nhpb).eq.0.0d0) then
2917 if (ibecarb(nhpb).eq.3) then
2918 jhpb(nhpb)=jhpb(nhpb)+nres
2919 else if (ibecarb(nhpb).eq.2) then
2920 ihpb(nhpb)=ihpb(nhpb)+nres
2921 else if (ibecarb(nhpb).eq.1) then
2922 ihpb(nhpb)=ihpb(nhpb)+nres
2923 jhpb(nhpb)=jhpb(nhpb)+nres
2925 if (dhpb(nhpb).eq.0.0d0)
2926 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2929 if (.not.out1file .or. me.eq.king)
2930 & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2931 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2933 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2934 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2937 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2938 C if (forcon(nhpb+1).gt.0.0d0) then
2940 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2943 if (restr_type.eq.4) then
2944 write (iout,*) "The BFAC array"
2946 write (iout,'(i5,f10.5)') i,bfac(i)
2949 if (itype(i).eq.ntyp1) cycle
2951 if (itype(j).eq.ntyp1) cycle
2952 if (itype(i).eq.10) then
2957 if (itype(j).eq.10) then
2967 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2970 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2971 ihpb(nhpb)=i+nres*ii
2972 jhpb(nhpb)=j+nres*jj
2973 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2975 if (.not.out1file .or. me.eq.king) then
2976 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2977 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2978 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2982 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2983 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2984 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2994 if (restr_type.eq.5) then
2995 restr_on_coord=.true.
2997 if (itype(i).eq.ntyp1) cycle
2998 bfac(i)=(scal_bfac/bfac(i))**2
3007 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
3008 & fordepthmax=fordepth(i)
3011 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
3016 c-------------------------------------------------------------------------------
3017 subroutine read_constr_homology
3019 include 'DIMENSIONS'
3023 include 'COMMON.SETUP'
3024 include 'COMMON.CONTROL'
3025 include 'COMMON.HOMOLOGY'
3026 include 'COMMON.CHAIN'
3027 include 'COMMON.IOUNITS'
3029 include 'COMMON.QRESTR'
3030 include 'COMMON.GEO'
3031 include 'COMMON.INTERACT'
3032 include 'COMMON.NAMES'
3034 c For new homol impl
3036 include 'COMMON.VAR'
3039 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
3041 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
3042 c & sigma_odl_temp(maxres,maxres,max_template)
3044 character*24 model_ki_dist, model_ki_angle
3045 character*500 controlcard
3046 integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
3047 & ik,iistart,nres_temp
3050 logical liiflag,lfirst
3053 c FP - Nov. 2014 Temporary specifications for new vars
3055 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
3057 double precision, dimension (max_template,maxres) :: rescore
3058 double precision, dimension (max_template,maxres) :: rescore2
3059 double precision, dimension (max_template,maxres) :: rescore3
3060 double precision distal
3061 character*24 tpl_k_rescore
3062 character*256 pdbfile
3063 c -----------------------------------------------------------------
3064 c Reading multiple PDB ref structures and calculation of retraints
3065 c not using pre-computed ones stored in files model_ki_{dist,angle}
3067 c -----------------------------------------------------------------
3070 c Alternative: reading from input
3071 call card_concat(controlcard)
3072 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
3073 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
3074 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
3075 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
3076 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
3077 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
3078 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3079 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3080 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3081 if(.not.read2sigma.and.start_from_model) then
3082 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3083 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3084 start_from_model=.false.
3085 iranconf=(indpdb.le.0)
3087 if(start_from_model .and. (me.eq.king .or. .not. out1file))
3088 & write(iout,*) 'START_FROM_MODELS is ON'
3089 c if(start_from_model .and. rest) then
3090 c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3091 c write(iout,*) 'START_FROM_MODELS is OFF'
3092 c write(iout,*) 'remove restart keyword from input'
3095 if (start_from_model) nmodel_start=constr_homology
3096 if (homol_nset.gt.1)then
3097 call card_concat(controlcard)
3098 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3099 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3100 write(iout,*) "iset homology_weight "
3102 write(iout,*) i,waga_homology(i)
3105 iset=mod(kolor,homol_nset)+1
3108 waga_homology(1)=1.0
3111 cd write (iout,*) "nnt",nnt," nct",nct
3118 c write(iout,*) 'nnt=',nnt,'nct=',nct
3121 do k=1,constr_homology
3134 if (read_homol_frag) then
3135 call read_klapaucjusz
3138 do k=1,constr_homology
3140 read(inp,'(a)') pdbfile
3141 pdbfiles_chomo(k)=pdbfile
3142 if(me.eq.king .or. .not. out1file)
3143 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3144 & pdbfile(:ilen(pdbfile))
3145 open(ipdbin,file=pdbfile,status='old',err=33)
3147 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3148 & pdbfile(:ilen(pdbfile))
3151 c print *,'Begin reading pdb data'
3153 c Files containing res sim or local scores (former containing sigmas)
3156 write(kic2,'(bz,i2.2)') k
3158 tpl_k_rescore="template"//kic2//".sco"
3162 if (read2sigma) then
3163 call readpdb_template(k)
3170 c Distance restraints
3173 C Copy the coordinates from reference coordinates (?)
3174 do i=1,2*nres_chomo(k)
3177 c write (iout,*) "c(",j,i,") =",c(j,i)
3181 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3183 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3184 open (ientin,file=tpl_k_rescore,status='old')
3185 if (nnt.gt.1) rescore(k,1)=0.0d0
3186 do irec=nnt,nct ! loop for reading res sim
3187 if (read2sigma) then
3188 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3189 & rescore3_tmp,idomain_tmp
3191 idomain(k,i_tmp)=idomain_tmp
3192 rescore(k,i_tmp)=rescore_tmp
3193 rescore2(k,i_tmp)=rescore2_tmp
3194 rescore3(k,i_tmp)=rescore3_tmp
3195 if (.not. out1file .or. me.eq.king)
3196 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3197 & i_tmp,rescore2_tmp,rescore_tmp,
3198 & rescore3_tmp,idomain_tmp
3201 read (ientin,*,end=1401) rescore_tmp
3203 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3204 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3205 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3210 if (waga_dist.ne.0.0d0) then
3218 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3219 c write (iout,*) k,i,j,distal,dist2_cut
3221 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3222 & .and. distal.le.dist2_cut ) then
3228 c write (iout,*) "k",k
3229 c write (iout,*) "i",i," j",j," constr_homology",
3234 if (read2sigma) then
3237 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3239 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3240 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3241 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3243 if (odl(k,ii).le.dist_cut) then
3244 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3247 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3248 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3250 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3251 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3255 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3258 l_homo(k,ii)=.false.
3264 c write (iout,*) "Distance restraints set"
3267 c Theta, dihedral and SC retraints
3269 if (waga_angle.gt.0.0d0) then
3270 c open (ientin,file=tpl_k_sigma_dih,status='old')
3271 c do irec=1,maxres-3 ! loop for reading sigma_dih
3272 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3273 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3274 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3275 c & sigma_dih(k,i+nnt-1)
3280 if (idomain(k,i).eq.0) then
3284 dih(k,i)=phiref(i) ! right?
3285 c read (ientin,*) sigma_dih(k,i) ! original variant
3286 c write (iout,*) "dih(",k,i,") =",dih(k,i)
3287 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3288 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3289 c & "rescore(",k,i-2,") =",rescore(k,i-2),
3290 c & "rescore(",k,i-3,") =",rescore(k,i-3)
3292 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3293 & rescore(k,i-2)+rescore(k,i-3))/4.0
3294 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3295 c write (iout,*) "Raw sigmas for dihedral angle restraints"
3296 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3297 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3298 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
3299 c Instead of res sim other local measure of b/b str reliability possible
3300 if (sigma_dih(k,i).ne.0)
3301 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3302 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3306 c write (iout,*) "Dihedral angle restraints set"
3309 if (waga_theta.gt.0.0d0) then
3310 c open (ientin,file=tpl_k_sigma_theta,status='old')
3311 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3312 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3313 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3314 c & sigma_theta(k,i+nnt-1)
3319 do i = nnt+2,nct ! right? without parallel.
3320 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3321 c do i=ithet_start,ithet_end ! with FG parallel.
3322 if (idomain(k,i).eq.0) then
3323 sigma_theta(k,i)=0.0
3326 thetatpl(k,i)=thetaref(i)
3327 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3328 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3329 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3330 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3331 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3332 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3333 & rescore(k,i-2))/3.0
3334 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3335 if (sigma_theta(k,i).ne.0)
3336 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3338 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3339 c rescore(k,i-2) ! right expression ?
3340 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3343 c write (iout,*) "Angle restraints set"
3346 if (waga_d.gt.0.0d0) then
3347 c open (ientin,file=tpl_k_sigma_d,status='old')
3348 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3349 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3350 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3351 c & sigma_d(k,i+nnt-1)
3355 do i = nnt,nct ! right? without parallel.
3356 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3357 c do i=loc_start,loc_end ! with FG parallel.
3358 if (itype(i).eq.10) cycle
3359 if (idomain(k,i).eq.0 ) then
3366 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3367 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3368 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3369 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3370 sigma_d(k,i)=rescore3(k,i) ! right expression ?
3371 if (sigma_d(k,i).ne.0)
3372 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3374 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3375 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3376 c read (ientin,*) sigma_d(k,i) ! 1st variant
3380 c write (iout,*) "SC restraints set"
3383 c remove distance restraints not used in any model from the list
3384 c shift data in all arrays
3386 c write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
3387 if (waga_dist.ne.0.0d0) then
3394 c if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3395 c & .and. distal.le.dist2_cut ) then
3396 c write (iout,*) "i",i," j",j," ii",ii
3398 if (ii_in_use(ii).eq.0.and.liiflag.or.
3399 & ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
3406 if(i10.eq.lim_odl) i10=i10+1
3408 ires_homo(iistart+ki)=ires_homo(ki+i01)
3409 jres_homo(iistart+ki)=jres_homo(ki+i01)
3410 ii_in_use(iistart+ki)=ii_in_use(ki+i01)
3411 do k=1,constr_homology
3412 odl(k,iistart+ki)=odl(k,ki+i01)
3413 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
3414 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
3417 iistart=iistart+i10-i01
3420 if (ii_in_use(ii).ne.0.and..not.liiflag) then
3428 c write (iout,*) "Removing distances completed"
3430 endif ! .not. klapaucjusz
3432 if (constr_homology.gt.0) call homology_partition
3433 c write (iout,*) "After homology_partition"
3435 if (constr_homology.gt.0) call init_int_table
3436 c write (iout,*) "After init_int_table"
3438 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3439 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3443 if (.not.out_template_restr) return
3444 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3445 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3446 write (iout,*) "Distance restraints from templates"
3448 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3449 & ii,ires_homo(ii),jres_homo(ii),
3450 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3451 & ki=1,constr_homology)
3453 write (iout,*) "Dihedral angle restraints from templates"
3455 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3456 & (rad2deg*dih(ki,i),
3457 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3459 write (iout,*) "Virtual-bond angle restraints from templates"
3461 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3462 & (rad2deg*thetatpl(ki,i),
3463 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3465 write (iout,*) "SC restraints from templates"
3467 write(iout,'(i5,100(4f8.2,4x))') i,
3468 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3469 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3472 c -----------------------------------------------------------------
3475 c----------------------------------------------------------------------
3477 subroutine flush(iu)
3482 subroutine flush(iu)
3487 c------------------------------------------------------------------------------
3488 subroutine copy_to_tmp(source)
3490 include "DIMENSIONS"
3491 include "COMMON.IOUNITS"
3492 character*(*) source
3493 character* 256 tmpfile
3497 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3498 inquire(file=tmpfile,exist=ex)
3500 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3501 & " to temporary directory..."
3502 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3503 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3507 c------------------------------------------------------------------------------
3508 subroutine move_from_tmp(source)
3510 include "DIMENSIONS"
3511 include "COMMON.IOUNITS"
3512 character*(*) source
3515 write (*,*) "Moving ",source(:ilen(source)),
3516 & " from temporary directory to working directory"
3517 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3518 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3521 c------------------------------------------------------------------------------
3522 subroutine random_init(seed)
3524 C Initialize random number generator
3527 include 'DIMENSIONS'
3530 logical OKRandom, prng_restart
3532 integer iseed_array(4)
3533 integer error_msg,ierr
3535 include 'COMMON.IOUNITS'
3536 include 'COMMON.TIME1'
3537 include 'COMMON.THREAD'
3538 include 'COMMON.SBRIDGE'
3539 include 'COMMON.CONTROL'
3540 include 'COMMON.MCM'
3541 include 'COMMON.MAP'
3542 include 'COMMON.HEADER'
3543 include 'COMMON.CSA'
3544 include 'COMMON.CHAIN'
3545 include 'COMMON.MUCA'
3547 include 'COMMON.FFIELD'
3548 include 'COMMON.SETUP'
3550 double precision seed,ran_number
3551 iseed=-dint(dabs(seed))
3552 if (iseed.eq.0) then
3553 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3554 & 'Random seed undefined. The program will stop.'
3555 write (*,'(/80(1h*)/20x,a/80(1h*))')
3556 & 'Random seed undefined. The program will stop.'
3558 call mpi_finalize(mpi_comm_world,ierr)
3560 stop 'Bad random seed.'
3563 if (fg_rank.eq.0) then
3567 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3568 OKRandom = prng_restart(me,iseed)
3571 tmp=65536.0d0**(4-i)
3572 iseed_array(i) = dint(seed/tmp)
3573 seed=seed-iseed_array(i)*tmp
3576 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3577 & (iseed_array(i),i=1,4)
3578 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3579 & (iseed_array(i),i=1,4)
3580 OKRandom = prng_restart(me,iseed_array)
3583 c r1 = prng_next(me)
3584 r1=ran_number(0.0D0,1.0D0)
3586 & write (iout,*) 'ran_num',r1
3587 if (r1.lt.0.0d0) OKRandom=.false.
3589 if (.not.OKRandom) then
3590 write (iout,*) 'PRNG IS NOT WORKING!!!'
3591 print *,'PRNG IS NOT WORKING!!!'
3594 call mpi_abort(mpi_comm_world,error_msg,ierr)
3597 write (iout,*) 'too many processors for parallel prng'
3598 write (*,*) 'too many processors for parallel prng'
3606 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3610 c----------------------------------------------------------------------
3611 subroutine read_klapaucjusz
3613 include 'DIMENSIONS'
3617 include 'COMMON.SETUP'
3618 include 'COMMON.CONTROL'
3619 include 'COMMON.HOMOLOGY'
3620 include 'COMMON.CHAIN'
3621 include 'COMMON.IOUNITS'
3623 include 'COMMON.GEO'
3624 include 'COMMON.INTERACT'
3625 include 'COMMON.NAMES'
3626 character*256 fragfile
3627 integer ninclust(maxclust),inclust(max_template,maxclust),
3628 & nresclust(maxclust),iresclust(maxres,maxclust),nclust
3631 character*24 model_ki_dist, model_ki_angle
3632 character*500 controlcard
3633 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp,
3634 & ik,ll,ii,kk,iistart,iishift,lim_xx
3635 double precision distal
3636 logical lprn /.true./
3643 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3644 double precision, dimension (max_template,maxres) :: rescore
3645 double precision, dimension (max_template,maxres) :: rescore2
3646 character*24 tpl_k_rescore
3647 character*256 pdbfile
3650 c For new homol impl
3652 include 'COMMON.VAR'
3654 call getenv("FRAGFILE",fragfile)
3655 open(ientin,file=fragfile,status="old",err=10)
3656 read(ientin,*) constr_homology,nclust
3657 nmodel_start=constr_homology
3663 do k=1,constr_homology
3664 read(ientin,'(a)') pdbfile
3665 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3666 & pdbfile(:ilen(pdbfile))
3667 pdbfiles_chomo(k)=pdbfile
3668 open(ipdbin,file=pdbfile,status='old',err=33)
3670 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3671 & pdbfile(:ilen(pdbfile))
3676 call readpdb_template(k)
3686 read(ientin,*) ninclust(i),nresclust(i)
3687 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3688 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3691 c Loop over clusters
3694 do ll = 1,ninclust(l)
3702 idomain(k,iresclust(i,l)+1) = 1
3704 idomain(k,iresclust(i,l)) = 1
3708 c Distance restraints
3711 C Copy the coordinates from reference coordinates (?)
3717 c write (iout,*) "c(",j,i,") =",c(j,i)
3720 call int_from_cart(.true.,.false.)
3721 call sc_loc_geom(.false.)
3723 thetaref(i)=theta(i)
3727 if (waga_dist.ne.0.0d0) then
3735 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3736 c write (iout,*) k,i,j,distal,dist2_cut
3738 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3739 & .and. distal.le.dist2_cut ) then
3745 c write (iout,*) "k",k
3746 c write (iout,*) "i",i," j",j," constr_homology",
3751 if (read2sigma) then
3754 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3756 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3757 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3758 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3760 if (odl(k,ii).le.dist_cut) then
3761 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3764 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3765 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3767 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3768 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3772 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3775 c l_homo(k,ii)=.false.
3782 c Theta, dihedral and SC retraints
3784 if (waga_angle.gt.0.0d0) then
3786 if (idomain(k,i).eq.0) then
3787 c sigma_dih(k,i)=0.0
3791 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3792 & rescore(k,i-2)+rescore(k,i-3))/4.0
3793 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3794 c & " sigma_dihed",sigma_dih(k,i)
3795 if (sigma_dih(k,i).ne.0)
3796 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3801 if (waga_theta.gt.0.0d0) then
3803 if (idomain(k,i).eq.0) then
3804 c sigma_theta(k,i)=0.0
3807 thetatpl(k,i)=thetaref(i)
3808 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3809 & rescore(k,i-2))/3.0
3810 if (sigma_theta(k,i).ne.0)
3811 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3815 if (waga_d.gt.0.0d0) then
3817 if (itype(i).eq.10) cycle
3818 if (idomain(k,i).eq.0 ) then
3825 sigma_d(k,i)=rescore(k,i)
3826 if (sigma_d(k,i).ne.0)
3827 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3828 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3834 c remove distance restraints not used in any model from the list
3835 c shift data in all arrays
3837 if (waga_dist.ne.0.0d0) then
3843 if (ii_in_use(ii).eq.0.and.liiflag) then
3847 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3848 & .not.liiflag.and.ii.eq.lim_odl) then
3849 if (ii.eq.lim_odl) then
3850 iishift=ii-iistart+1
3855 do ki=iistart,lim_odl-iishift
3856 ires_homo(ki)=ires_homo(ki+iishift)
3857 jres_homo(ki)=jres_homo(ki+iishift)
3858 ii_in_use(ki)=ii_in_use(ki+iishift)
3859 do k=1,constr_homology
3860 odl(k,ki)=odl(k,ki+iishift)
3861 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3862 l_homo(k,ki)=l_homo(k,ki+iishift)
3866 lim_odl=lim_odl-iishift
3873 10 stop "Error in fragment file"