8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors
534 write (iout,*) 'FTORS factor =',ftors
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
556 if ( secstruc(i) .eq. 'H') then
557 ! Helix restraints for this residue
560 phi0(ii) = 45.0D0*deg2rad
561 drange(ii)= 5.0D0*deg2rad
562 phibound(1,i) = phi0(ii)-drange(ii)
563 phibound(2,i) = phi0(ii)+drange(ii)
564 else if (secstruc(i) .eq. 'E') then
565 ! strand restraints for this residue
568 phi0(ii) = 180.0D0*deg2rad
569 drange(ii)= 5.0D0*deg2rad
570 phibound(1,i) = phi0(ii)-drange(ii)
571 phibound(2,i) = phi0(ii)+drange(ii)
573 ! no restraints for this residue
574 ndih_nconstr=ndih_nconstr+1
575 idih_nconstr(ndih_nconstr)=i
579 ! deallocate(secstruc)
582 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
583 ! deallocate(secstruc)
586 write(iout,'(A20)')'Error reading FTORS'
587 ! deallocate(secstruc)
589 end subroutine secstrp2dihc
590 !-----------------------------------------------------------------------------
591 subroutine read_secstr_pred(jin,jout,errors)
593 ! implicit real*8 (a-h,o-z)
594 ! INCLUDE 'DIMENSIONS'
595 ! include 'COMMON.IOUNITS'
596 ! include 'COMMON.CHAIN'
597 !el character(len=1),dimension(nres) :: secstruc !(maxres)
598 !el COMMON/SECONDARYS/secstruc
600 character(len=80) :: line,line1 !,ucase
601 logical :: errflag,errors,blankline
604 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
607 read (jin,'(a)') line
608 write (jout,'(2a)') '> ',line(1:78)
610 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
611 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
612 ! end-groups in the input file "*.spred"
615 do while (index(line1,'$END').eq.0)
616 ! Override commented lines.
619 do while (.not.blankline)
621 call mykey(line,line1,ipos,blankline,errflag)
622 if (errflag) write (jout,'(2a)') &
623 'Error when reading sequence in line: ',line
624 errors=errors .or. errflag
625 if (.not. blankline .and. .not. errflag) then
628 !el if (iseq.le.maxres) then
629 if (line1(1:1).eq.'-' ) then
630 secstruc(iseq)=line1(1:1)
631 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
632 ( ucase(line1(1:1)).eq.'H' ) ) then
633 secstruc(iseq)=ucase(line1(1:1))
636 write (jout,1010) line1(1:1), iseq
641 !el write (jout,1000) iseq,maxres
644 do while (ipos1.le.iend)
649 !el if (iseq.le.maxres) then
650 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
651 secstruc(iseq)=line1(ipos1-1:ipos1-1)
652 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
653 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
654 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
657 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
662 !el write (jout,1000) iseq,maxres
669 read (jin,'(a)') line
670 write (jout,'(2a)') '> ',line(1:78)
674 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
676 !d check whether the found length of the chain is correct.
677 length_of_chain=iseq-1
678 if (length_of_chain .ne. nres) then
680 write (jout,'(a,i4,a,i4,a)') &
681 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
682 ,length_of_chain,') does not match with the number of residues (' &
687 1000 format('Error - the number of residues (',i4,&
688 ') has exceeded maximum (',i4,').')
689 1010 format ('Error - unrecognized secondary structure label',a4,&
692 end subroutine read_secstr_pred
694 !-----------------------------------------------------------------------------
696 !-----------------------------------------------------------------------------
701 use control_data, only:maxterm !,maxtor
705 use control, only: getenv_loc
707 ! Read the parameters of the probability distributions of the virtual-bond
708 ! valence angles and the side chains and energy parameters.
710 ! Important! Energy-term weights ARE NOT read here; they are read from the
711 ! main input file instead, because NO defaults have yet been set for these
714 ! implicit real*8 (a-h,o-z)
715 ! include 'DIMENSIONS'
720 ! include 'COMMON.IOUNITS'
721 ! include 'COMMON.CHAIN'
722 ! include 'COMMON.INTERACT'
723 ! include 'COMMON.GEO'
724 ! include 'COMMON.LOCAL'
725 ! include 'COMMON.TORSION'
726 ! include 'COMMON.SCCOR'
727 ! include 'COMMON.SCROT'
728 ! include 'COMMON.FFIELD'
729 ! include 'COMMON.NAMES'
730 ! include 'COMMON.SBRIDGE'
731 ! include 'COMMON.MD'
732 ! include 'COMMON.SETUP'
733 character(len=1) :: t1,t2,t3
734 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
735 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
736 logical :: lprint,LaTeX
737 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
738 real(kind=8),dimension(13) :: b
739 character(len=3) :: lancuch !,ucase
741 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
742 integer :: maxinter,junk,kk,ii
743 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
744 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
745 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 integer :: ichir1,ichir2
748 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
749 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
750 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
752 ! For printing parameters after they are read set the following in the UNRES
755 ! setenv PRINT_PARM YES
757 ! To print parameters in LaTeX format rather than as ASCII tables:
761 call getenv_loc("PRINT_PARM",lancuch)
762 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
763 call getenv_loc("LATEX",lancuch)
764 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
766 dwa16=2.0d0**(1.0d0/6.0d0)
768 ! Assign virtual-bond length
771 vblinv2=vblinv*vblinv
773 ! Read the virtual-bond parameters, masses, and moments of inertia
774 ! and Stokes' radii of the peptide group and side chains
776 allocate(dsc(ntyp1)) !(ntyp1)
777 allocate(dsc_inv(ntyp1)) !(ntyp1)
778 allocate(nbondterm(ntyp)) !(ntyp)
779 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
780 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
781 allocate(msc(ntyp+1)) !(ntyp+1)
782 allocate(isc(ntyp+1)) !(ntyp+1)
783 allocate(restok(ntyp+1)) !(ntyp+1)
784 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
790 read (ibond,*) vbldp0,akp,mp,ip,pstok
793 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
794 dsc(i) = vbldsc0(1,i)
798 dsc_inv(i)=1.0D0/dsc(i)
802 read (ibond,*) junk,vbldpDUM,vbldp0,akp,rjunk,mp,ip,pstok
804 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
805 j=1,nbondterm(i)),msc(i),isc(i),restok(i)
806 dsc(i) = vbldsc0(1,i)
810 dsc_inv(i)=1.0D0/dsc(i)
815 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
816 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
818 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
820 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),&
821 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
823 write (iout,'(13x,3f10.5)') &
824 vbldsc0(j,i),aksc(j,i),abond0(j,i)
828 !----------------------------------------------------
829 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
830 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
831 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
832 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
833 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
834 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
847 ! Read the parameters of the probability distribution/energy expression
848 ! of the virtual-bond valence angles theta
851 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
852 (bthet(j,i,1,1),j=1,2)
853 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
854 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
855 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
859 athet(1,i,1,-1)=athet(1,i,1,1)
860 athet(2,i,1,-1)=athet(2,i,1,1)
861 bthet(1,i,1,-1)=-bthet(1,i,1,1)
862 bthet(2,i,1,-1)=-bthet(2,i,1,1)
863 athet(1,i,-1,1)=-athet(1,i,1,1)
864 athet(2,i,-1,1)=-athet(2,i,1,1)
865 bthet(1,i,-1,1)=bthet(1,i,1,1)
866 bthet(2,i,-1,1)=bthet(2,i,1,1)
870 athet(1,i,-1,-1)=athet(1,-i,1,1)
871 athet(2,i,-1,-1)=-athet(2,-i,1,1)
872 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
873 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
874 athet(1,i,-1,1)=athet(1,-i,1,1)
875 athet(2,i,-1,1)=-athet(2,-i,1,1)
876 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
877 bthet(2,i,-1,1)=bthet(2,-i,1,1)
878 athet(1,i,1,-1)=-athet(1,-i,1,1)
879 athet(2,i,1,-1)=athet(2,-i,1,1)
880 bthet(1,i,1,-1)=bthet(1,-i,1,1)
881 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
886 polthet(j,i)=polthet(j,-i)
889 gthet(j,i)=gthet(j,-i)
897 'Parameters of the virtual-bond valence angles:'
898 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
899 ' ATHETA0 ',' A1 ',' A2 ',&
902 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
903 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
905 write (iout,'(/a/9x,5a/79(1h-))') &
906 'Parameters of the expression for sigma(theta_c):',&
907 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
908 ' ALPH3 ',' SIGMA0C '
910 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
911 (polthet(j,i),j=0,3),sigc0(i)
913 write (iout,'(/a/9x,5a/79(1h-))') &
914 'Parameters of the second gaussian:',&
915 ' THETA0 ',' SIGMA0 ',' G1 ',&
918 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),&
919 sig0(i),(gthet(j,i),j=1,3)
923 'Parameters of the virtual-bond valence angles:'
924 write (iout,'(/a/9x,5a/79(1h-))') &
925 'Coefficients of expansion',&
926 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
927 ' b1*10^1 ',' b2*10^1 '
929 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
930 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
931 (10*bthet(j,i,1,1),j=1,2)
933 write (iout,'(/a/9x,5a/79(1h-))') &
934 'Parameters of the expression for sigma(theta_c):',&
935 ' alpha0 ',' alph1 ',' alph2 ',&
936 ' alhp3 ',' sigma0c '
938 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
939 (polthet(j,i),j=0,3),sigc0(i)
941 write (iout,'(/a/9x,5a/79(1h-))') &
942 'Parameters of the second gaussian:',&
943 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
946 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
947 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
953 ! Read the parameters of Utheta determined from ab initio surfaces
954 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
956 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
957 ntheterm3,nsingle,ndouble
958 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
960 !----------------------------------------------------
961 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
962 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
963 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
964 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
965 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
966 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
967 !(maxtheterm,-maxthetyp1:maxthetyp1,&
968 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
969 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
970 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
971 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
972 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
973 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
974 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
975 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
976 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
977 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
978 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
979 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
980 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
981 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
982 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
983 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
984 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
986 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
988 ithetyp(i)=-ithetyp(-i)
991 aa0thet(:,:,:,:)=0.0d0
992 aathet(:,:,:,:,:)=0.0d0
993 bbthet(:,:,:,:,:,:)=0.0d0
994 ccthet(:,:,:,:,:,:)=0.0d0
995 ddthet(:,:,:,:,:,:)=0.0d0
996 eethet(:,:,:,:,:,:)=0.0d0
997 ffthet(:,:,:,:,:,:,:)=0.0d0
998 ggthet(:,:,:,:,:,:,:)=0.0d0
1000 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1002 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1003 ! VAR:1=non-glicyne non-proline 2=proline
1004 ! VAR:negative values for D-aminoacid
1006 do j=-nthetyp,nthetyp
1007 do k=-nthetyp,nthetyp
1008 read (ithep,'(6a)',end=111,err=111) res1
1009 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1010 ! VAR: aa0thet is variable describing the average value of Foureir
1011 ! VAR: expansion series
1012 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1013 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1014 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1015 read (ithep,*,end=111,err=111) &
1016 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1017 read (ithep,*,end=111,err=111) &
1018 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1019 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1020 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1021 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1023 read (ithep,*,end=111,err=111) &
1024 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1025 ffthet(lll,llll,ll,i,j,k,iblock),&
1026 ggthet(llll,lll,ll,i,j,k,iblock),&
1027 ggthet(lll,llll,ll,i,j,k,iblock),&
1028 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1033 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1034 ! coefficients of theta-and-gamma-dependent terms are zero.
1035 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1036 ! RECOMENTDED AFTER VERSION 3.3)
1040 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1041 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1043 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1044 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1047 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1049 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1052 ! AND COMMENT THE LOOPS BELOW
1056 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1057 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1059 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1060 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1063 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1065 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1070 ! Substitution for D aminoacids from symmetry.
1073 do j=-nthetyp,nthetyp
1074 do k=-nthetyp,nthetyp
1075 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1077 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1081 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1082 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1083 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1084 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1090 ffthet(llll,lll,ll,i,j,k,iblock)= &
1091 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1092 ffthet(lll,llll,ll,i,j,k,iblock)= &
1093 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1094 ggthet(llll,lll,ll,i,j,k,iblock)= &
1095 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1096 ggthet(lll,llll,ll,i,j,k,iblock)= &
1097 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1106 ! Control printout of the coefficients of virtual-bond-angle potentials
1109 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1114 write (iout,'(//4a)') &
1115 'Type ',onelett(i),onelett(j),onelett(k)
1116 write (iout,'(//a,10x,a)') " l","a[l]"
1117 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1118 write (iout,'(i2,1pe15.5)') &
1119 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1121 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1122 "b",l,"c",l,"d",l,"e",l
1124 write (iout,'(i2,4(1pe15.5))') m,&
1125 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1126 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1130 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1131 "f+",l,"f-",l,"g+",l,"g-",l
1134 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1135 ffthet(n,m,l,i,j,k,iblock),&
1136 ffthet(m,n,l,i,j,k,iblock),&
1137 ggthet(n,m,l,i,j,k,iblock),&
1138 ggthet(m,n,l,i,j,k,iblock)
1148 write (2,*) "Start reading THETA_PDB",ithep_pdb
1150 ! write (2,*) 'i=',i
1151 read (ithep_pdb,*,err=111,end=111) &
1152 a0thet(i),(athet(j,i,1,1),j=1,2),&
1153 (bthet(j,i,1,1),j=1,2)
1154 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1155 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1156 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1157 sigc0(i)=sigc0(i)**2
1160 athet(1,i,1,-1)=athet(1,i,1,1)
1161 athet(2,i,1,-1)=athet(2,i,1,1)
1162 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1163 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1164 athet(1,i,-1,1)=-athet(1,i,1,1)
1165 athet(2,i,-1,1)=-athet(2,i,1,1)
1166 bthet(1,i,-1,1)=bthet(1,i,1,1)
1167 bthet(2,i,-1,1)=bthet(2,i,1,1)
1170 a0thet(i)=a0thet(-i)
1171 athet(1,i,-1,-1)=athet(1,-i,1,1)
1172 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1173 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1174 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1175 athet(1,i,-1,1)=athet(1,-i,1,1)
1176 athet(2,i,-1,1)=-athet(2,-i,1,1)
1177 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1178 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1179 athet(1,i,1,-1)=-athet(1,-i,1,1)
1180 athet(2,i,1,-1)=athet(2,-i,1,1)
1181 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1182 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1183 theta0(i)=theta0(-i)
1187 polthet(j,i)=polthet(j,-i)
1190 gthet(j,i)=gthet(j,-i)
1193 write (2,*) "End reading THETA_PDB"
1198 !-------------------------------------------
1199 allocate(nlob(ntyp1)) !(ntyp1)
1200 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1201 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1202 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1209 gaussc(:,:,:,:)=0.0D0
1213 ! Read the parameters of the probability distribution/energy expression
1214 ! of the side chains.
1217 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1221 dsc_inv(i)=1.0D0/dsc(i)
1232 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1233 ((blower(k,l,1),l=1,k),k=1,3)
1234 censc(1,1,-i)=censc(1,1,i)
1235 censc(2,1,-i)=censc(2,1,i)
1236 censc(3,1,-i)=-censc(3,1,i)
1238 read (irotam,*,end=112,err=112) bsc(j,i)
1239 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1240 ((blower(k,l,j),l=1,k),k=1,3)
1241 censc(1,j,-i)=censc(1,j,i)
1242 censc(2,j,-i)=censc(2,j,i)
1243 censc(3,j,-i)=-censc(3,j,i)
1244 ! BSC is amplitude of Gaussian
1251 akl=akl+blower(k,m,j)*blower(l,m,j)
1255 if (((k.eq.3).and.(l.ne.3)) &
1256 .or.((l.eq.3).and.(k.ne.3))) then
1257 gaussc(k,l,j,-i)=-akl
1258 gaussc(l,k,j,-i)=-akl
1260 gaussc(k,l,j,-i)=akl
1261 gaussc(l,k,j,-i)=akl
1270 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1273 if (nlobi.gt.0) then
1275 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
1276 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1277 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1278 'log h',(bsc(j,i),j=1,nlobi)
1279 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1280 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1282 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1283 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1286 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1287 write (iout,'(a,f10.4,4(16x,f10.4))') &
1288 'Center ',(bsc(j,i),j=1,nlobi)
1289 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1298 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1299 ! added by Urszula Kozlowska 07/11/2007
1301 !el Maximum number of SC local term fitting function coefficiants
1302 !el integer,parameter :: maxsccoef=65
1304 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1307 read (irotam,*,end=112,err=112)
1309 read (irotam,*,end=112,err=112)
1312 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1317 ! Read the parameters of the probability distribution/energy expression
1318 ! of the side chains.
1320 write (2,*) "Start reading ROTAM_PDB"
1322 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1326 dsc_inv(i)=1.0D0/dsc(i)
1337 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1338 ((blower(k,l,1),l=1,k),k=1,3)
1340 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1341 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1342 ((blower(k,l,j),l=1,k),k=1,3)
1349 akl=akl+blower(k,m,j)*blower(l,m,j)
1359 write (2,*) "End reading ROTAM_PDB"
1365 ! Read torsional parameters in old format
1367 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1369 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1370 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1371 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1373 !el from energy module--------
1374 allocate(v1(nterm_old,ntortyp,ntortyp))
1375 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1376 !el---------------------------
1381 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1387 write (iout,'(/a/)') 'Torsional constants:'
1390 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1391 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1397 ! Read torsional parameters
1399 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1400 read (itorp,*,end=113,err=113) ntortyp
1401 !el from energy module---------
1402 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1403 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1405 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1406 allocate(vlor2(maxlor,ntortyp,ntortyp))
1407 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1408 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1410 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1411 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1412 !el---------------------------
1415 !el---------------------------
1417 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1419 itortyp(i)=-itortyp(-i)
1421 ! itortyp(ntyp1)=ntortyp
1422 ! itortyp(-ntyp1)=-ntortyp
1424 write (iout,*) 'ntortyp',ntortyp
1426 do j=-ntortyp+1,ntortyp-1
1427 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1429 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1430 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1433 do k=1,nterm(i,j,iblock)
1434 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1436 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1437 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1438 v0ij=v0ij+si*v1(k,i,j,iblock)
1440 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1441 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1442 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1444 do k=1,nlor(i,j,iblock)
1445 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1446 vlor2(k,i,j),vlor3(k,i,j)
1447 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1450 v0(-i,-j,iblock)=v0ij
1456 write (iout,'(/a/)') 'Torsional constants:'
1458 do i=-ntortyp,ntortyp
1459 do j=-ntortyp,ntortyp
1460 write (iout,*) 'ityp',i,' jtyp',j
1461 write (iout,*) 'Fourier constants'
1462 do k=1,nterm(i,j,iblock)
1463 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1466 write (iout,*) 'Lorenz constants'
1467 do k=1,nlor(i,j,iblock)
1468 write (iout,'(3(1pe15.5))') &
1469 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1475 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1477 ! 6/23/01 Read parameters for double torsionals
1479 !el from energy module------------
1480 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1481 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1482 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1483 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1484 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1485 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1486 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1487 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1488 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1489 !---------------------------------
1493 do j=-ntortyp+1,ntortyp-1
1494 do k=-ntortyp+1,ntortyp-1
1495 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1496 ! write (iout,*) "OK onelett",
1499 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1500 .or. t3.ne.toronelet(k)) then
1501 write (iout,*) "Error in double torsional parameter file",&
1504 call MPI_Finalize(Ierror)
1506 stop "Error in double torsional parameter file"
1508 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1509 ntermd_2(i,j,k,iblock)
1510 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1511 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1512 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1513 ntermd_1(i,j,k,iblock))
1514 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1515 ntermd_1(i,j,k,iblock))
1516 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1517 ntermd_1(i,j,k,iblock))
1518 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1519 ntermd_1(i,j,k,iblock))
1520 ! Martix of D parameters for one dimesional foureir series
1521 do l=1,ntermd_1(i,j,k,iblock)
1522 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1523 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1524 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1525 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1526 ! write(iout,*) "whcodze" ,
1527 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1529 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1530 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1531 v2s(m,l,i,j,k,iblock),&
1532 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1533 ! Martix of D parameters for two dimesional fourier series
1534 do l=1,ntermd_2(i,j,k,iblock)
1536 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1537 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1538 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1539 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1548 write (iout,*) 'Constants for double torsionals'
1551 do j=-ntortyp+1,ntortyp-1
1552 do k=-ntortyp+1,ntortyp-1
1553 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1554 ' nsingle',ntermd_1(i,j,k,iblock),&
1555 ' ndouble',ntermd_2(i,j,k,iblock)
1557 write (iout,*) 'Single angles:'
1558 do l=1,ntermd_1(i,j,k,iblock)
1559 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1560 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1561 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1562 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1565 write (iout,*) 'Pairs of angles:'
1566 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1567 do l=1,ntermd_2(i,j,k,iblock)
1568 write (iout,'(i5,20f10.5)') &
1569 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1572 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1573 do l=1,ntermd_2(i,j,k,iblock)
1574 write (iout,'(i5,20f10.5)') &
1575 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1576 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1585 ! Read of Side-chain backbone correlation parameters
1586 ! Modified 11 May 2012 by Adasko
1589 read (isccor,*,end=119,err=119) nsccortyp
1591 !el from module energy-------------
1592 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1593 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1594 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1595 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1596 !-----------------------------------
1598 !el from module energy-------------
1599 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1601 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1603 isccortyp(i)=-isccortyp(-i)
1605 iscprol=isccortyp(20)
1606 ! write (iout,*) 'ntortyp',ntortyp
1608 !c maxinter is maximum interaction sites
1609 !el from module energy---------
1610 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1611 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1612 -nsccortyp:nsccortyp))
1613 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1614 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1615 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1616 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1617 !-----------------------------------
1621 read (isccor,*,end=119,err=119) &
1622 nterm_sccor(i,j),nlor_sccor(i,j)
1628 nterm_sccor(-i,j)=nterm_sccor(i,j)
1629 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1630 nterm_sccor(i,-j)=nterm_sccor(i,j)
1631 do k=1,nterm_sccor(i,j)
1632 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1634 if (j.eq.iscprol) then
1635 if (i.eq.isccortyp(10)) then
1636 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1637 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1639 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1640 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1641 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1642 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1643 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1644 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1645 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1646 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1649 if (i.eq.isccortyp(10)) then
1650 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1651 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1653 if (j.eq.isccortyp(10)) then
1654 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1655 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1657 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1658 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1659 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1660 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1661 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1662 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1666 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1667 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1668 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1669 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1672 do k=1,nlor_sccor(i,j)
1673 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1674 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1675 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1676 (1+vlor3sccor(k,i,j)**2)
1678 v0sccor(l,i,j)=v0ijsccor
1679 v0sccor(l,-i,j)=v0ijsccor1
1680 v0sccor(l,i,-j)=v0ijsccor2
1681 v0sccor(l,-i,-j)=v0ijsccor3
1687 !el from module energy-------------
1688 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1690 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1691 ! write (iout,*) 'ntortyp',ntortyp
1693 !c maxinter is maximum interaction sites
1694 !el from module energy---------
1695 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1696 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1697 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1698 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1699 !-----------------------------------
1703 read (isccor,*,end=119,err=119) &
1704 nterm_sccor(i,j),nlor_sccor(i,j)
1708 do k=1,nterm_sccor(i,j)
1709 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1711 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1714 do k=1,nlor_sccor(i,j)
1715 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1716 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1717 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1718 (1+vlor3sccor(k,i,j)**2)
1720 v0sccor(l,i,j)=v0ijsccor !el ,iblock
1728 write (iout,'(/a/)') 'Torsional constants:'
1731 write (iout,*) 'ityp',i,' jtyp',j
1732 write (iout,*) 'Fourier constants'
1733 do k=1,nterm_sccor(i,j)
1734 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1736 write (iout,*) 'Lorenz constants'
1737 do k=1,nlor_sccor(i,j)
1738 write (iout,'(3(1pe15.5))') &
1739 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1746 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1747 ! interaction energy of the Gly, Ala, and Pro prototypes.
1752 write (iout,*) "Coefficients of the cumulants"
1754 read (ifourier,*) nloctyp
1755 !write(iout,*) "nloctyp",nloctyp
1756 !el from module energy-------
1757 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1758 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1759 allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1760 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1761 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1762 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1763 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1764 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1768 !--------------------------------
1771 read (ifourier,*,end=115,err=115)
1772 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
1774 write (iout,*) 'Type',i
1775 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
1785 B1tilde(1,-i) =-b(3)
1787 ! b1tilde(1,i)=0.0d0
1788 ! b1tilde(2,i)=0.0d0
1813 Ctilde(1,2,-i)=-b(9)
1817 ! Ctilde(1,1,i)=0.0d0
1818 ! Ctilde(1,2,i)=0.0d0
1819 ! Ctilde(2,1,i)=0.0d0
1820 ! Ctilde(2,2,i)=0.0d0
1838 Dtilde(1,2,-i)=-b(8)
1842 ! Dtilde(1,1,i)=0.0d0
1843 ! Dtilde(1,2,i)=0.0d0
1844 ! Dtilde(2,1,i)=0.0d0
1845 ! Dtilde(2,2,i)=0.0d0
1846 EE(1,1,i)= b(10)+b(11)
1847 EE(2,2,i)=-b(10)+b(11)
1848 EE(2,1,i)= b(12)-b(13)
1849 EE(1,2,i)= b(12)+b(13)
1850 EE(1,1,-i)= b(10)+b(11)
1851 EE(2,2,-i)=-b(10)+b(11)
1852 EE(2,1,-i)=-b(12)+b(13)
1853 EE(1,2,-i)=-b(12)-b(13)
1859 ! ee(2,1,i)=ee(1,2,i)
1863 write (iout,*) 'Type',i
1865 write(iout,*) B1(1,i),B1(2,i)
1867 write(iout,*) B2(1,i),B2(2,i)
1870 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1874 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1878 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1883 ! Read electrostatic-interaction parameters
1888 write (iout,'(/a)') 'Electrostatic interaction constants:'
1889 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
1890 'IT','JT','APP','BPP','AEL6','AEL3'
1892 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1893 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1894 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1895 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1900 app (i,j)=epp(i,j)*rri*rri
1901 bpp (i,j)=-2.0D0*epp(i,j)*rri
1902 ael6(i,j)=elpp6(i,j)*4.2D0**6
1903 ael3(i,j)=elpp3(i,j)*4.2D0**3
1905 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
1911 ! Read side-chain interaction parameters.
1913 !el from module energy - COMMON.INTERACT-------
1914 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
1915 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
1916 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
1917 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
1918 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
1927 !--------------------------------
1929 read (isidep,*,end=117,err=117) ipot,expon
1930 if (ipot.lt.1 .or. ipot.gt.5) then
1931 write (iout,'(2a)') 'Error while reading SC interaction',&
1932 'potential file - unknown potential type.'
1934 call MPI_Finalize(Ierror)
1940 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
1941 ', exponents are ',expon,2*expon
1942 ! goto (10,20,30,30,40) ipot
1944 !----------------------- LJ potential ---------------------------------
1946 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1947 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1948 (sigma0(i),i=1,ntyp)
1950 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1951 write (iout,'(a/)') 'The epsilon array:'
1952 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1953 write (iout,'(/a)') 'One-body parameters:'
1954 write (iout,'(a,4x,a)') 'residue','sigma'
1955 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1958 !----------------------- LJK potential --------------------------------
1960 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1961 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1962 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1964 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1965 write (iout,'(a/)') 'The epsilon array:'
1966 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1967 write (iout,'(/a)') 'One-body parameters:'
1968 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1969 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
1973 !---------------------- GB or BP potential -----------------------------
1977 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1979 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1980 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1981 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1982 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1983 ! For the GB potential convert sigma'**2 into chi'
1986 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
1990 write (iout,'(/a/)') 'Parameters of the BP potential:'
1991 write (iout,'(a/)') 'The epsilon array:'
1992 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1993 write (iout,'(/a)') 'One-body parameters:'
1994 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
1996 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
1997 chip(i),alp(i),i=1,ntyp)
2000 !--------------------- GBV potential -----------------------------------
2002 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2003 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2004 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2005 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2007 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2008 write (iout,'(a/)') 'The epsilon array:'
2009 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2010 write (iout,'(/a)') 'One-body parameters:'
2011 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2012 's||/s_|_^2',' chip ',' alph '
2013 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
2014 sigii(i),chip(i),alp(i),i=1,ntyp)
2017 write(iout,*)"Wrong ipot"
2022 !-----------------------------------------------------------------------
2023 ! Calculate the "working" parameters of SC interactions.
2025 !el from module energy - COMMON.INTERACT-------
2026 allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2027 allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2034 !--------------------------------
2043 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2044 sigma(j,i)=sigma(i,j)
2045 rs0(i,j)=dwa16*sigma(i,j)
2049 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2050 'Working parameters of the SC interactions:',&
2051 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2056 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2065 sigeps=dsign(1.0D0,epsij)
2067 aa(i,j)=epsij*rrij*rrij
2068 bb(i,j)=-sigeps*epsij*rrij
2072 sigt1sq=sigma0(i)**2
2073 sigt2sq=sigma0(j)**2
2076 ratsig1=sigt2sq/sigt1sq
2077 ratsig2=1.0D0/ratsig1
2078 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2079 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2080 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2084 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2085 sigmaii(i,j)=rsum_max
2086 sigmaii(j,i)=rsum_max
2088 ! sigmaii(i,j)=r0(i,j)
2089 ! sigmaii(j,i)=r0(i,j)
2091 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2092 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2093 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2094 augm(i,j)=epsij*r_augm**(2*expon)
2095 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2102 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2103 restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),&
2104 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2109 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2114 ! Define the SC-p interaction constants (hard-coded; old style)
2117 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2119 ! aad(i,1)=0.3D0*4.0D0**12
2120 ! Following line for constants currently implemented
2121 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2122 aad(i,1)=1.5D0*4.0D0**12
2123 ! aad(i,1)=0.17D0*5.6D0**12
2125 ! "Soft" SC-p repulsion
2127 ! Following line for constants currently implemented
2128 ! aad(i,1)=0.3D0*4.0D0**6
2129 ! "Hard" SC-p repulsion
2130 bad(i,1)=3.0D0*4.0D0**6
2131 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2140 ! 8/9/01 Read the SC-p interaction constants from file
2143 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2146 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2147 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2148 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2149 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2153 write (iout,*) "Parameters of SC-p interactions:"
2155 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2156 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2162 ! Define the constants of the disulfide bridge
2166 ! Old arbitrary potential - commented out.
2171 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2172 ! energy surface of diethyl disulfide.
2173 ! A. Liwo and U. Kozlowska, 11/24/03
2190 write (iout,'(/a)') "Disulfide bridge parameters:"
2191 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2192 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2193 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2194 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2198 111 write (iout,*) "Error reading bending energy parameters."
2200 112 write (iout,*) "Error reading rotamer energy parameters."
2202 113 write (iout,*) "Error reading torsional energy parameters."
2204 114 write (iout,*) "Error reading double torsional energy parameters."
2206 115 write (iout,*) &
2207 "Error reading cumulant (multibody energy) parameters."
2209 116 write (iout,*) "Error reading electrostatic energy parameters."
2211 117 write (iout,*) "Error reading side chain interaction parameters."
2213 118 write (iout,*) "Error reading SCp interaction parameters."
2215 119 write (iout,*) "Error reading SCCOR parameters"
2218 call MPI_Finalize(Ierror)
2222 end subroutine parmread
2224 !-----------------------------------------------------------------------------
2226 !-----------------------------------------------------------------------------
2227 subroutine printmat(ldim,m,n,iout,key,a)
2230 character(len=3),dimension(n) :: key
2231 real(kind=8),dimension(ldim,n) :: a
2233 integer :: i,j,k,m,iout,nlim
2237 write (iout,1000) (key(k),k=i,nlim)
2239 1000 format (/5x,8(6x,a3))
2240 1020 format (/80(1h-)/)
2242 write (iout,1010) key(j),(a(j,k),k=i,nlim)
2245 1010 format (a3,2x,8(f9.4))
2247 end subroutine printmat
2248 !-----------------------------------------------------------------------------
2250 !-----------------------------------------------------------------------------
2252 ! Read the PDB file and convert the peptide geometry into virtual-chain
2255 use energy_data, only: itype
2259 use control, only: rescode
2260 ! implicit real*8 (a-h,o-z)
2261 ! include 'DIMENSIONS'
2262 ! include 'COMMON.LOCAL'
2263 ! include 'COMMON.VAR'
2264 ! include 'COMMON.CHAIN'
2265 ! include 'COMMON.INTERACT'
2266 ! include 'COMMON.IOUNITS'
2267 ! include 'COMMON.GEO'
2268 ! include 'COMMON.NAMES'
2269 ! include 'COMMON.CONTROL'
2270 ! include 'COMMON.DISTFIT'
2271 ! include 'COMMON.SETUP'
2272 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift!,ity!,&
2274 logical :: lprn=.true.,fail
2275 real(kind=8),dimension(3) :: e1,e2,e3
2276 real(kind=8) :: dcj,efree_temp
2277 character(len=3) :: seq,res
2278 character(len=5) :: atom
2279 character(len=80) :: card
2280 real(kind=8),dimension(3,20) :: sccor
2281 integer :: kkk,lll,icha,kupa !rescode,
2284 integer,dimension(2,maxres/3) :: hfrag_alloc
2285 integer,dimension(4,maxres/3) :: bfrag_alloc
2286 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2292 ! write (2,*) "UNRES_PDB",unres_pdb
2300 !-----------------------------
2301 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2302 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2305 read (ipdbin,'(a80)',end=10) card
2306 ! write (iout,'(a)') card
2307 if (card(:5).eq.'HELIX') then
2310 read(card(22:25),*) hfrag(1,nhfrag)
2311 read(card(34:37),*) hfrag(2,nhfrag)
2313 if (card(:5).eq.'SHEET') then
2316 read(card(24:26),*) bfrag(1,nbfrag)
2317 read(card(35:37),*) bfrag(2,nbfrag)
2318 !rc----------------------------------------
2319 !rc to be corrected !!!
2320 bfrag(3,nbfrag)=bfrag(1,nbfrag)
2321 bfrag(4,nbfrag)=bfrag(2,nbfrag)
2322 !rc----------------------------------------
2324 if (card(:3).eq.'END') then
2326 else if (card(:3).eq.'TER') then
2330 itype(ires_old)=ntyp1
2331 itype(ires_old-1)=ntyp1
2333 ! write (iout,*) "Chain ended",ires,ishift,ires_old
2336 dc(j,ires)=sccor(j,iii)
2339 call sccenter(ires,iii,sccor)
2345 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2346 ! Fish out the ATOM cards.
2347 if (index(card(1:4),'ATOM').gt.0) then
2348 read (card(12:16),*) atom
2349 ! write (iout,*) "! ",atom," !",ires
2350 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
2351 read (card(23:26),*) ires
2352 read (card(18:20),'(a3)') res
2353 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
2354 ! & " ires_old",ires_old
2355 ! write (iout,*) "ishift",ishift," ishift1",ishift1
2356 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2357 if (ires-ishift+ishift1.ne.ires_old) then
2358 ! Calculate the CM of the preceding residue.
2359 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2361 ! write (iout,*) "Calculating sidechain center iii",iii
2364 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2367 call sccenter(ires_old,iii,sccor)
2371 ! Start new residue.
2372 if (res.eq.'Cl-' .or. res.eq.'Na+') then
2375 else if (ibeg.eq.1) then
2376 write (iout,*) "BEG ires",ires
2378 if (res.ne.'GLY' .and. res.ne. 'ACE') then
2382 ires=ires-ishift+ishift1
2384 ! write (iout,*) "ishift",ishift," ires",ires,&
2385 ! " ires_old",ires_old
2387 else if (ibeg.eq.2) then
2389 ishift=-ires_old+ires-1 !!!!!
2390 ishift1=ishift1-1 !!!!!
2391 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2392 ires=ires-ishift+ishift1
2396 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2397 ires=ires-ishift+ishift1
2400 if (res.eq.'ACE' .or. res.eq.'NHE') then
2403 itype(ires)=rescode(ires,res,0)
2406 ires=ires-ishift+ishift1
2408 ! write (iout,*) "ires_old",ires_old," ires",ires
2409 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2412 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
2413 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2414 res.eq.'NHE'.and.atom(:2).eq.'HN') then
2415 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2416 ! write (iout,*) "backbone ",atom
2418 write (iout,'(2i3,2x,a,3f8.3)') &
2419 ires,itype(ires),res,(c(j,ires),j=1,3)
2423 sccor(j,iii)=c(j,ires)
2425 ! write (*,*) card(23:27),ires,itype(ires)
2426 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2427 atom.ne.'N' .and. atom.ne.'C' .and. &
2428 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2429 atom.ne.'OXT' .and. atom(:2).ne.'3H') then
2430 ! write (iout,*) "sidechain ",atom
2432 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2436 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
2437 if (ires.eq.0) return
2438 ! Calculate dummy residue coordinates inside the "chain" of a multichain
2442 ! write (iout,*) i,itype(i)
2443 ! if (itype(i).eq.ntyp1) then
2444 ! write (iout,*) "dummy",i,itype(i)
2446 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
2447 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
2451 if (itype(i).eq.ntyp1) then
2452 if (itype(i+1).eq.ntyp1) then
2453 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
2454 ! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
2455 ! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
2457 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2458 ! print *,i,'tu dochodze'
2459 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
2467 c(j,i)=c(j,i-1)-1.9d0*e2(j)
2471 dcj=(c(j,i-2)-c(j,i-3))/2.0
2472 if (dcj.eq.0) dcj=1.23591524223
2477 else !itype(i+1).eq.ntyp1
2479 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2480 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
2487 c(j,i)=c(j,i+1)-1.9d0*e2(j)
2491 dcj=(c(j,i+3)-c(j,i+2))/2.0
2492 if (dcj.eq.0) dcj=1.23591524223
2497 endif !itype(i+1).eq.ntyp1
2498 endif !itype.eq.ntyp1
2501 ! Calculate the CM of the last side chain.
2505 dc(j,ires)=sccor(j,iii)
2508 call sccenter(ires,iii,sccor)
2514 if (itype(nres).ne.10) then
2518 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2519 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
2526 c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
2530 dcj=c(j,nres-2)-c(j,nres-3)
2531 c(j,nres)=c(j,nres-1)+dcj
2532 c(j,2*nres)=c(j,nres)
2536 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
2538 if (nres.ne.nres0) then
2539 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
2541 stop "Error nres value in WHAM input"
2544 !---------------------------------
2545 !el reallocate tables
2548 ! hfrag_alloc(j,i)=hfrag(j,i)
2551 ! bfrag_alloc(j,i)=bfrag(j,i)
2557 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
2558 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
2559 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
2560 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
2564 ! hfrag(j,i)=hfrag_alloc(j,i)
2569 ! bfrag(j,i)=bfrag_alloc(j,i)
2572 !el end reallocate tables
2573 !---------------------------------
2581 c(j,2*nres)=c(j,nres)
2583 if (itype(1).eq.ntyp1) then
2587 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2588 call refsys(2,3,4,e1,e2,e3,fail)
2595 c(j,1)=c(j,2)-3.8d0*e2(j)
2605 ! Copy the coordinates to reference coordinates
2611 ! Calculate internal coordinates.
2613 write (iout,'(/a)') &
2614 "Cartesian coordinates of the reference structure"
2615 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
2616 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
2618 write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
2619 restyp(itype(ires)),ires,(c(j,ires),j=1,3),&
2620 (c(j,ires+nres),j=1,3)
2623 ! znamy już nres wiec mozna alokowac tablice
2624 ! Calculate internal coordinates.
2625 if(me.eq.king.or..not.out1file)then
2626 write (iout,'(a)') &
2627 "Backbone and SC coordinates as read from the PDB"
2629 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
2630 ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),&
2631 (c(j,nres+ires),j=1,3)
2635 if(.not.allocated(vbld)) then
2636 allocate(vbld(2*nres))
2641 if(.not.allocated(vbld_inv)) then
2642 allocate(vbld_inv(2*nres))
2648 if(.not.allocated(theta)) then
2649 allocate(theta(nres+2))
2653 if(.not.allocated(phi)) allocate(phi(nres+2))
2654 if(.not.allocated(alph)) allocate(alph(nres+2))
2655 if(.not.allocated(omeg)) allocate(omeg(nres+2))
2656 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
2657 if(.not.allocated(phiref)) allocate(phiref(nres+2))
2658 if(.not.allocated(costtab)) allocate(costtab(nres))
2659 if(.not.allocated(sinttab)) allocate(sinttab(nres))
2660 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
2661 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
2662 if(.not.allocated(xxref)) allocate(xxref(nres))
2663 if(.not.allocated(yyref)) allocate(yyref(nres))
2664 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
2665 if(.not.allocated(dc_norm)) then
2666 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
2667 allocate(dc_norm(3,0:2*nres+2))
2671 call int_from_cart(.true.,.false.)
2672 call sc_loc_geom(.false.)
2674 thetaref(i)=theta(i)
2684 dc(j,i)=c(j,i+1)-c(j,i)
2685 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2690 dc(j,i+nres)=c(j,i+nres)-c(j,i)
2691 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2693 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
2697 ! Copy the coordinates to reference coordinates
2698 ! Splits to single chain if occurs
2702 ! cref(j,i,cou)=c(j,i)
2706 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
2707 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
2708 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
2709 !-----------------------------
2715 !c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2717 if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
2720 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2725 cref(j,i,cou)=c(j,i)
2726 cref(j,i+nres,cou)=c(j,i+nres)
2728 chain_rep(j,lll,kkk)=c(j,i)
2729 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
2733 write (iout,*) chain_length
2734 if (chain_length.eq.0) chain_length=nres
2736 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
2737 chain_rep(j,chain_length+nres,symetr) &
2738 =chain_rep(j,chain_length+nres,1)
2741 ! write (iout,*) "spraw lancuchy",chain_length,symetr
2743 ! do kkk=1,chain_length
2744 ! write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
2748 ! makes copy of chains
2749 write (iout,*) "symetr", symetr
2751 if (symetr.gt.1) then
2758 write(iout,*) (tabperm(i,kkk),kkk=1,4)
2764 ! write (iout,*) i,icha
2765 do lll=1,chain_length
2767 if (cou.le.nres) then
2769 kupa=mod(lll,chain_length)
2770 iprzes=(kkk-1)*chain_length+lll
2771 if (kupa.eq.0) kupa=chain_length
2772 ! write (iout,*) "kupa", kupa
2773 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
2774 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
2781 !-koniec robienia kopii
2784 write (iout,*) "nowa struktura", nperm
2786 write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),&
2788 cref(3,i,kkk),cref(1,nres+i,kkk),&
2789 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
2791 100 format (//' alpha-carbon coordinates ',&
2792 ' centroid coordinates'/ &
2793 ' ', 6X,'X',11X,'Y',11X,'Z', &
2794 10X,'X',11X,'Y',11X,'Z')
2795 110 format (a,'(',i3,')',6f12.5)
2801 bfrag(i,j)=bfrag(i,j)-ishift
2807 hfrag(i,j)=hfrag(i,j)-ishift
2813 end subroutine readpdb
2814 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
2815 !-----------------------------------------------------------------------------
2817 !-----------------------------------------------------------------------------
2818 subroutine read_control
2832 use random, only: random_init
2833 ! implicit real*8 (a-h,o-z)
2834 ! include 'DIMENSIONS'
2836 use prng, only:prng_restart
2838 logical :: OKRandom!, prng_restart
2841 ! include 'COMMON.IOUNITS'
2842 ! include 'COMMON.TIME1'
2843 ! include 'COMMON.THREAD'
2844 ! include 'COMMON.SBRIDGE'
2845 ! include 'COMMON.CONTROL'
2846 ! include 'COMMON.MCM'
2847 ! include 'COMMON.MAP'
2848 ! include 'COMMON.HEADER'
2849 ! include 'COMMON.CSA'
2850 ! include 'COMMON.CHAIN'
2851 ! include 'COMMON.MUCA'
2852 ! include 'COMMON.MD'
2853 ! include 'COMMON.FFIELD'
2854 ! include 'COMMON.INTERACT'
2855 ! include 'COMMON.SETUP'
2856 !el integer :: KDIAG,ICORFL,IXDR
2857 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
2858 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
2859 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
2860 ! character(len=80) :: ucase
2861 character(len=640) :: controlcard
2863 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
2869 read (INP,'(a)') titel
2870 call card_concat(controlcard,.true.)
2871 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
2872 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
2873 call reada(controlcard,'SEED',seed,0.0D0)
2874 call random_init(seed)
2875 ! Set up the time limit (caution! The time must be input in minutes!)
2876 read_cart=index(controlcard,'READ_CART').gt.0
2877 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
2878 call readi(controlcard,'SYM',symetr,1)
2879 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
2880 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
2881 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
2882 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
2883 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
2884 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
2885 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
2886 call reada(controlcard,'DRMS',drms,0.1D0)
2887 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2888 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
2889 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
2890 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
2891 write (iout,'(a,f10.1)')'DRMS = ',drms
2892 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
2893 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
2895 call readi(controlcard,'NZ_START',nz_start,0)
2896 call readi(controlcard,'NZ_END',nz_end,0)
2897 call readi(controlcard,'IZ_SC',iz_sc,0)
2898 timlim=60.0D0*timlim
2899 safety = 60.0d0*safety
2902 call reada(controlcard,"T_BATH",t_bath,300.0d0)
2903 !C Varibles set size of box
2904 call reada(controlcard,'BOXX',boxxsize,100.0d0)
2905 call reada(controlcard,'BOXY',boxysize,100.0d0)
2906 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
2907 ! CUTOFFF ON ELECTROSTATICS
2908 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
2909 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
2911 !C-------------------------
2912 minim=(index(controlcard,'MINIMIZE').gt.0)
2913 dccart=(index(controlcard,'CART').gt.0)
2914 overlapsc=(index(controlcard,'OVERLAP').gt.0)
2915 overlapsc=.not.overlapsc
2916 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
2917 searchsc=.not.searchsc
2918 sideadd=(index(controlcard,'SIDEADD').gt.0)
2919 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
2920 outpdb=(index(controlcard,'PDBOUT').gt.0)
2921 outmol2=(index(controlcard,'MOL2OUT').gt.0)
2922 pdbref=(index(controlcard,'PDBREF').gt.0)
2923 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
2924 indpdb=index(controlcard,'PDBSTART')
2925 extconf=(index(controlcard,'EXTCONF').gt.0)
2926 call readi(controlcard,'IPRINT',iprint,0)
2927 call readi(controlcard,'MAXGEN',maxgen,10000)
2928 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
2929 call readi(controlcard,"KDIAG",kdiag,0)
2930 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
2931 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
2932 write (iout,*) "RESCALE_MODE",rescale_mode
2933 split_ene=index(controlcard,'SPLIT_ENE').gt.0
2934 if (index(controlcard,'REGULAR').gt.0.0D0) then
2935 call reada(controlcard,'WEIDIS',weidis,0.1D0)
2939 if (index(controlcard,'CHECKGRAD').gt.0) then
2941 if (index(controlcard,'CART').gt.0) then
2943 elseif (index(controlcard,'CARINT').gt.0) then
2948 elseif (index(controlcard,'THREAD').gt.0) then
2950 call readi(controlcard,'THREAD',nthread,0)
2951 if (nthread.gt.0) then
2952 call reada(controlcard,'WEIDIS',weidis,0.1D0)
2955 write (iout,'(a)')'A number has to follow the THREAD keyword.'
2956 stop 'Error termination in Read_Control.'
2958 else if (index(controlcard,'MCMA').gt.0) then
2960 else if (index(controlcard,'MCEE').gt.0) then
2962 else if (index(controlcard,'MULTCONF').gt.0) then
2964 else if (index(controlcard,'MAP').gt.0) then
2966 call readi(controlcard,'MAP',nmap,0)
2967 else if (index(controlcard,'CSA').gt.0) then
2969 !rc else if (index(controlcard,'ZSCORE').gt.0) then
2971 !rc ZSCORE is rm from UNRES, modecalc=9 is available
2974 !fcm else if (index(controlcard,'MCMF').gt.0) then
2976 else if (index(controlcard,'SOFTREG').gt.0) then
2978 else if (index(controlcard,'CHECK_BOND').gt.0) then
2980 else if (index(controlcard,'TEST').gt.0) then
2982 else if (index(controlcard,'MD').gt.0) then
2984 else if (index(controlcard,'RE ').gt.0) then
2988 lmuca=index(controlcard,'MUCA').gt.0
2989 call readi(controlcard,'MUCADYN',mucadyn,0)
2990 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
2991 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
2993 write (iout,*) 'MUCADYN=',mucadyn
2994 write (iout,*) 'MUCASMOOTH=',muca_smooth
2997 iscode=index(controlcard,'ONE_LETTER')
2998 indphi=index(controlcard,'PHI')
2999 indback=index(controlcard,'BACK')
3000 iranconf=index(controlcard,'RAND_CONF')
3001 i2ndstr=index(controlcard,'USE_SEC_PRED')
3002 gradout=index(controlcard,'GRADOUT').gt.0
3003 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3004 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3005 if (me.eq.king .or. .not.out1file ) &
3006 write (iout,*) "DISTCHAINMAX",distchainmax
3008 if(me.eq.king.or..not.out1file) &
3009 write (iout,'(2a)') diagmeth(kdiag),&
3010 ' routine used to diagonalize matrices.'
3012 end subroutine read_control
3013 !-----------------------------------------------------------------------------
3014 subroutine read_REMDpar
3016 ! Read REMD settings
3023 use control_data, only:out1file
3024 ! implicit real*8 (a-h,o-z)
3025 ! include 'DIMENSIONS'
3026 ! include 'COMMON.IOUNITS'
3027 ! include 'COMMON.TIME1'
3028 ! include 'COMMON.MD'
3031 !el include 'COMMON.LANGEVIN'
3033 !el include 'COMMON.LANGEVIN.lang0'
3035 ! include 'COMMON.INTERACT'
3036 ! include 'COMMON.NAMES'
3037 ! include 'COMMON.GEO'
3038 ! include 'COMMON.REMD'
3039 ! include 'COMMON.CONTROL'
3040 ! include 'COMMON.SETUP'
3041 ! character(len=80) :: ucase
3042 character(len=320) :: controlcard
3043 character(len=3200) :: controlcard1
3044 integer :: iremd_m_total
3047 ! real(kind=8) :: var,ene
3049 if(me.eq.king.or..not.out1file) &
3050 write (iout,*) "REMD setup"
3052 call card_concat(controlcard,.true.)
3053 call readi(controlcard,"NREP",nrep,3)
3054 call readi(controlcard,"NSTEX",nstex,1000)
3055 call reada(controlcard,"RETMIN",retmin,10.0d0)
3056 call reada(controlcard,"RETMAX",retmax,1000.0d0)
3057 mremdsync=(index(controlcard,'SYNC').gt.0)
3058 call readi(controlcard,"NSYN",i_sync_step,100)
3059 restart1file=(index(controlcard,'REST1FILE').gt.0)
3060 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3061 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3062 if(max_cache_traj_use.gt.max_cache_traj) &
3063 max_cache_traj_use=max_cache_traj
3064 if(me.eq.king.or..not.out1file) then
3065 !d if (traj1file) then
3066 !rc caching is in testing - NTWX is not ignored
3067 !d write (iout,*) "NTWX value is ignored"
3068 !d write (iout,*) " trajectory is stored to one file by master"
3069 !d write (iout,*) " before exchange at NSTEX intervals"
3071 write (iout,*) "NREP= ",nrep
3072 write (iout,*) "NSTEX= ",nstex
3073 write (iout,*) "SYNC= ",mremdsync
3074 write (iout,*) "NSYN= ",i_sync_step
3075 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3078 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3079 if (index(controlcard,'TLIST').gt.0) then
3081 call card_concat(controlcard1,.true.)
3082 read(controlcard1,*) (remd_t(i),i=1,nrep)
3083 if(me.eq.king.or..not.out1file) &
3084 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3087 if (index(controlcard,'MLIST').gt.0) then
3089 call card_concat(controlcard1,.true.)
3090 read(controlcard1,*) (remd_m(i),i=1,nrep)
3091 if(me.eq.king.or..not.out1file) then
3092 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3095 iremd_m_total=iremd_m_total+remd_m(i)
3097 write (iout,*) 'Total number of replicas ',iremd_m_total
3100 if(me.eq.king.or..not.out1file) &
3101 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3103 end subroutine read_REMDpar
3104 !-----------------------------------------------------------------------------
3105 subroutine read_MDpar
3109 use control_data, only: r_cut,rlamb,out1file
3111 use geometry_data, only: pi
3113 ! implicit real*8 (a-h,o-z)
3114 ! include 'DIMENSIONS'
3115 ! include 'COMMON.IOUNITS'
3116 ! include 'COMMON.TIME1'
3117 ! include 'COMMON.MD'
3120 !el include 'COMMON.LANGEVIN'
3122 !el include 'COMMON.LANGEVIN.lang0'
3124 ! include 'COMMON.INTERACT'
3125 ! include 'COMMON.NAMES'
3126 ! include 'COMMON.GEO'
3127 ! include 'COMMON.SETUP'
3128 ! include 'COMMON.CONTROL'
3129 ! include 'COMMON.SPLITELE'
3130 ! character(len=80) :: ucase
3131 character(len=320) :: controlcard
3136 call card_concat(controlcard,.true.)
3137 call readi(controlcard,"NSTEP",n_timestep,1000000)
3138 call readi(controlcard,"NTWE",ntwe,100)
3139 call readi(controlcard,"NTWX",ntwx,1000)
3140 call reada(controlcard,"DT",d_time,1.0d-1)
3141 call reada(controlcard,"DVMAX",dvmax,2.0d1)
3142 call reada(controlcard,"DAMAX",damax,1.0d1)
3143 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3144 call readi(controlcard,"LANG",lang,0)
3145 RESPA = index(controlcard,"RESPA") .gt. 0
3146 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3147 ntime_split0=ntime_split
3148 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3149 ntime_split0=ntime_split
3150 call reada(controlcard,"R_CUT",r_cut,2.0d0)
3151 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3152 rest = index(controlcard,"REST").gt.0
3153 tbf = index(controlcard,"TBF").gt.0
3154 usampl = index(controlcard,"USAMPL").gt.0
3155 mdpdb = index(controlcard,"MDPDB").gt.0
3156 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3157 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3158 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3159 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3160 if (count_reset_moment.eq.0) count_reset_moment=1000000000
3161 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3162 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3163 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3164 if (count_reset_vel.eq.0) count_reset_vel=1000000000
3165 large = index(controlcard,"LARGE").gt.0
3166 print_compon = index(controlcard,"PRINT_COMPON").gt.0
3167 rattle = index(controlcard,"RATTLE").gt.0
3168 ! if performing umbrella sampling, fragments constrained are read from the fragment file
3174 if(me.eq.king.or..not.out1file) then
3176 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3178 write (iout,'(a)') "The units are:"
3179 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3180 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3181 " acceleration: angstrom/(48.9 fs)**2"
3182 write (iout,'(a)') "energy: kcal/mol, temperature: K"
3184 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3185 write (iout,'(a60,f10.5,a)') &
3186 "Initial time step of numerical integration:",d_time,&
3188 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3190 write (iout,'(2a,i4,a)') &
3191 "A-MTS algorithm used; initial time step for fast-varying",&
3192 " short-range forces split into",ntime_split," steps."
3193 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3194 r_cut," lambda",rlamb
3196 write (iout,'(2a,f10.5)') &
3197 "Maximum acceleration threshold to reduce the time step",&
3198 "/increase split number:",damax
3199 write (iout,'(2a,f10.5)') &
3200 "Maximum predicted energy drift to reduce the timestep",&
3201 "/increase split number:",edriftmax
3202 write (iout,'(a60,f10.5)') &
3203 "Maximum velocity threshold to reduce velocities:",dvmax
3204 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3205 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3206 if (rattle) write (iout,'(a60)') &
3207 "Rattle algorithm used to constrain the virtual bonds"
3211 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3212 call reada(controlcard,"RWAT",rwat,1.4d0)
3213 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3214 surfarea=index(controlcard,"SURFAREA").gt.0
3215 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3216 if(me.eq.king.or..not.out1file)then
3217 write (iout,'(/a,$)') "Langevin dynamics calculation"
3219 write (iout,'(a/)') &
3220 " with direct integration of Langevin equations"
3221 else if (lang.eq.2) then
3222 write (iout,'(a/)') " with TINKER stochasic MD integrator"
3223 else if (lang.eq.3) then
3224 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3225 else if (lang.eq.4) then
3226 write (iout,'(a/)') " in overdamped mode"
3228 write (iout,'(//a,i5)') &
3229 "=========== ERROR: Unknown Langevin dynamics mode:",lang
3232 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3233 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
3234 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
3235 write (iout,'(a60,f10.5)') &
3236 "Scaling factor of the friction forces:",scal_fric
3237 if (surfarea) write (iout,'(2a,i10,a)') &
3238 "Friction coefficients will be scaled by solvent-accessible",&
3239 " surface area every",reset_fricmat," steps."
3241 ! Calculate friction coefficients and bounds of stochastic forces
3242 eta=6*pi*cPoise*etawat
3243 if(me.eq.king.or..not.out1file) &
3244 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
3246 gamp=scal_fric*(pstok+rwat)*eta
3247 stdfp=dsqrt(2*Rb*t_bath/d_time)
3248 allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1)
3250 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
3251 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
3253 if(me.eq.king.or..not.out1file)then
3254 write (iout,'(/2a/)') &
3255 "Radii of site types and friction coefficients and std's of",&
3256 " stochastic forces of fully exposed sites"
3257 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
3259 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),&
3260 gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
3264 if(me.eq.king.or..not.out1file)then
3265 write (iout,'(a)') "Berendsen bath calculation"
3266 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3267 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
3269 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
3270 count_reset_moment," steps"
3272 write (iout,'(a,i10,a)') &
3273 "Velocities will be reset at random every",count_reset_vel,&
3277 if(me.eq.king.or..not.out1file) &
3278 write (iout,'(a31)') "Microcanonical mode calculation"
3280 if(me.eq.king.or..not.out1file)then
3281 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
3283 write(iout,*) "MD running with constraints."
3284 write(iout,*) "Equilibration time ", eq_time, " mtus."
3285 write(iout,*) "Constraining ", nfrag," fragments."
3286 write(iout,*) "Length of each fragment, weight and q0:"
3288 write (iout,*) "Set of restraints #",iset
3290 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
3291 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
3293 write(iout,*) "constraints between ", npair, "fragments."
3294 write(iout,*) "constraint pairs, weights and q0:"
3296 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
3297 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
3299 write(iout,*) "angle constraints within ", nfrag_back,&
3300 "backbone fragments."
3301 write(iout,*) "fragment, weights:"
3303 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
3304 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
3305 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
3308 iset=mod(kolor,nset)+1
3311 if(me.eq.king.or..not.out1file) &
3312 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
3314 end subroutine read_MDpar
3315 !-----------------------------------------------------------------------------
3319 ! implicit real*8 (a-h,o-z)
3320 ! include 'DIMENSIONS'
3321 ! include 'COMMON.MAP'
3322 ! include 'COMMON.IOUNITS'
3323 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
3324 character(len=80) :: mapcard !,ucase
3327 ! real(kind=8) :: var,ene
3330 read (inp,'(a)') mapcard
3331 mapcard=ucase(mapcard)
3332 if (index(mapcard,'PHI').gt.0) then
3334 else if (index(mapcard,'THE').gt.0) then
3336 else if (index(mapcard,'ALP').gt.0) then
3338 else if (index(mapcard,'OME').gt.0) then
3341 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
3342 stop 'Error - illegal variable spec in MAP card.'
3344 call readi (mapcard,'RES1',res1(imap),0)
3345 call readi (mapcard,'RES2',res2(imap),0)
3346 if (res1(imap).eq.0) then
3347 res1(imap)=res2(imap)
3348 else if (res2(imap).eq.0) then
3349 res2(imap)=res1(imap)
3351 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
3352 write (iout,'(a)') &
3353 'Error - illegal definition of variable group in MAP.'
3354 stop 'Error - illegal definition of variable group in MAP.'
3356 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
3357 call reada(mapcard,'TO',ang_to(imap),0.0D0)
3358 call readi(mapcard,'NSTEP',nstep(imap),0)
3359 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
3360 write (iout,'(a)') &
3361 'Illegal boundary and/or step size specification in MAP.'
3362 stop 'Illegal boundary and/or step size specification in MAP.'
3366 end subroutine map_read
3367 !-----------------------------------------------------------------------------
3370 use control_data, only: vdisulf
3372 ! implicit real*8 (a-h,o-z)
3373 ! include 'DIMENSIONS'
3374 ! include 'COMMON.IOUNITS'
3375 ! include 'COMMON.GEO'
3376 ! include 'COMMON.CSA'
3377 ! include 'COMMON.BANK'
3378 ! include 'COMMON.CONTROL'
3379 ! character(len=80) :: ucase
3380 character(len=620) :: mcmcard
3382 ! integer :: ntf,ik,iw_pdb
3383 ! real(kind=8) :: var,ene
3385 call card_concat(mcmcard,.true.)
3387 call readi(mcmcard,'NCONF',nconf,50)
3388 call readi(mcmcard,'NADD',nadd,0)
3389 call readi(mcmcard,'JSTART',jstart,1)
3390 call readi(mcmcard,'JEND',jend,1)
3391 call readi(mcmcard,'NSTMAX',nstmax,500000)
3392 call readi(mcmcard,'N0',n0,1)
3393 call readi(mcmcard,'N1',n1,6)
3394 call readi(mcmcard,'N2',n2,4)
3395 call readi(mcmcard,'N3',n3,0)
3396 call readi(mcmcard,'N4',n4,0)
3397 call readi(mcmcard,'N5',n5,0)
3398 call readi(mcmcard,'N6',n6,10)
3399 call readi(mcmcard,'N7',n7,0)
3400 call readi(mcmcard,'N8',n8,0)
3401 call readi(mcmcard,'N9',n9,0)
3402 call readi(mcmcard,'N14',n14,0)
3403 call readi(mcmcard,'N15',n15,0)
3404 call readi(mcmcard,'N16',n16,0)
3405 call readi(mcmcard,'N17',n17,0)
3406 call readi(mcmcard,'N18',n18,0)
3408 vdisulf=(index(mcmcard,'DYNSS').gt.0)
3410 call readi(mcmcard,'NDIFF',ndiff,2)
3411 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
3412 call readi(mcmcard,'IS1',is1,1)
3413 call readi(mcmcard,'IS2',is2,8)
3414 call readi(mcmcard,'NRAN0',nran0,4)
3415 call readi(mcmcard,'NRAN1',nran1,2)
3416 call readi(mcmcard,'IRR',irr,1)
3417 call readi(mcmcard,'NSEED',nseed,20)
3418 call readi(mcmcard,'NTOTAL',ntotal,10000)
3419 call reada(mcmcard,'CUT1',cut1,2.0d0)
3420 call reada(mcmcard,'CUT2',cut2,5.0d0)
3421 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
3422 call readi(mcmcard,'ICMAX',icmax,3)
3423 call readi(mcmcard,'IRESTART',irestart,0)
3424 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
3427 call reada(mcmcard,'DELE',dele,20.0d0)
3428 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
3429 call readi(mcmcard,'IREF',iref,0)
3430 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
3431 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
3432 call readi(mcmcard,'NCONF_IN',nconf_in,0)
3433 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
3434 write (iout,*) "NCONF_IN",nconf_in
3436 end subroutine csaread
3437 !-----------------------------------------------------------------------------
3441 use control_data, only: MaxMoveType
3444 ! implicit real*8 (a-h,o-z)
3445 ! include 'DIMENSIONS'
3446 ! include 'COMMON.MCM'
3447 ! include 'COMMON.MCE'
3448 ! include 'COMMON.IOUNITS'
3449 ! character(len=80) :: ucase
3450 character(len=320) :: mcmcard
3453 ! real(kind=8) :: var,ene
3455 call card_concat(mcmcard,.true.)
3456 call readi(mcmcard,'MAXACC',maxacc,100)
3457 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
3458 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
3459 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
3460 call readi(mcmcard,'MAXREPM',maxrepm,200)
3461 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
3462 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
3463 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
3464 call reada(mcmcard,'E_UP',e_up,5.0D0)
3465 call reada(mcmcard,'DELTE',delte,0.1D0)
3466 call readi(mcmcard,'NSWEEP',nsweep,5)
3467 call readi(mcmcard,'NSTEPH',nsteph,0)
3468 call readi(mcmcard,'NSTEPC',nstepc,0)
3469 call reada(mcmcard,'TMIN',tmin,298.0D0)
3470 call reada(mcmcard,'TMAX',tmax,298.0D0)
3471 call readi(mcmcard,'NWINDOW',nwindow,0)
3472 call readi(mcmcard,'PRINT_MC',print_mc,0)
3473 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
3474 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
3475 ent_read=(index(mcmcard,'ENT_READ').gt.0)
3476 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
3477 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
3478 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
3479 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
3480 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
3481 if (nwindow.gt.0) then
3482 allocate(winstart(nwindow)) !!el (maxres)
3483 allocate(winend(nwindow)) !!el
3484 allocate(winlen(nwindow)) !!el
3485 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
3487 winlen(i)=winend(i)-winstart(i)+1
3490 if (tmax.lt.tmin) tmax=tmin
3491 if (tmax.eq.tmin) then
3495 if (nstepc.gt.0 .and. nsteph.gt.0) then
3496 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
3497 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
3499 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
3500 ! Probabilities of different move types
3501 sumpro_type(0)=0.0D0
3502 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
3503 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
3504 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
3505 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
3506 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
3507 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
3508 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
3510 print *,'i',i,' sumprotype',sumpro_type(i)
3511 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
3512 print *,'i',i,' sumprotype',sumpro_type(i)
3515 end subroutine mcmread
3516 !-----------------------------------------------------------------------------
3517 subroutine read_minim
3520 ! implicit real*8 (a-h,o-z)
3521 ! include 'DIMENSIONS'
3522 ! include 'COMMON.MINIM'
3523 ! include 'COMMON.IOUNITS'
3524 ! character(len=80) :: ucase
3525 character(len=320) :: minimcard
3527 ! integer :: ntf,ik,iw_pdb
3528 ! real(kind=8) :: var,ene
3530 call card_concat(minimcard,.true.)
3531 call readi(minimcard,'MAXMIN',maxmin,2000)
3532 call readi(minimcard,'MAXFUN',maxfun,5000)
3533 call readi(minimcard,'MINMIN',minmin,maxmin)
3534 call readi(minimcard,'MINFUN',minfun,maxmin)
3535 call reada(minimcard,'TOLF',tolf,1.0D-2)
3536 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
3537 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
3538 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
3539 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
3540 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
3541 'Options in energy minimization:'
3542 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
3543 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
3544 'MinMin:',MinMin,' MinFun:',MinFun,&
3545 ' TolF:',TolF,' RTolF:',RTolF
3547 end subroutine read_minim
3548 !-----------------------------------------------------------------------------
3549 subroutine openunits
3551 use MD_data, only: usampl
3554 use control_data, only:out1file
3555 use control, only: getenv_loc
3556 ! implicit real*8 (a-h,o-z)
3557 ! include 'DIMENSIONS'
3560 character(len=16) :: form,nodename
3561 integer :: nodelen,ierror,npos
3563 ! include 'COMMON.SETUP'
3564 ! include 'COMMON.IOUNITS'
3565 ! include 'COMMON.MD'
3566 ! include 'COMMON.CONTROL'
3567 integer :: lenpre,lenpot,lentmp !,ilen
3569 character(len=3) :: out1file_text !,ucase
3570 character(len=3) :: ll
3573 ! integer :: ntf,ik,iw_pdb
3574 ! real(kind=8) :: var,ene
3576 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
3577 call getenv_loc("PREFIX",prefix)
3579 call getenv_loc("POT",pot)
3580 call getenv_loc("DIRTMP",tmpdir)
3581 call getenv_loc("CURDIR",curdir)
3582 call getenv_loc("OUT1FILE",out1file_text)
3583 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
3584 out1file_text=ucase(out1file_text)
3585 if (out1file_text(1:1).eq."Y") then
3588 out1file=fg_rank.gt.0
3593 if (lentmp.gt.0) then
3594 write (*,'(80(1h!))')
3595 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
3596 write (*,'(80(1h!))')
3597 write (*,*)"All output files will be on node /tmp directory."
3599 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
3600 if (me.eq.king) then
3601 write (*,*) "The master node is ",nodename
3602 else if (fg_rank.eq.0) then
3603 write (*,*) "I am the CG slave node ",nodename
3605 write (*,*) "I am the FG slave node ",nodename
3608 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
3609 lenpre = lentmp+lenpre+1
3611 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
3612 ! Get the names and open the input files
3613 #if defined(WINIFL) || defined(WINPGI)
3614 open(1,file=pref_orig(:ilen(pref_orig))// &
3615 '.inp',status='old',readonly,shared)
3616 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3617 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3618 ! Get parameter filenames and open the parameter files.
3619 call getenv_loc('BONDPAR',bondname)
3620 open (ibond,file=bondname,status='old',readonly,shared)
3621 call getenv_loc('THETPAR',thetname)
3622 open (ithep,file=thetname,status='old',readonly,shared)
3623 call getenv_loc('ROTPAR',rotname)
3624 open (irotam,file=rotname,status='old',readonly,shared)
3625 call getenv_loc('TORPAR',torname)
3626 open (itorp,file=torname,status='old',readonly,shared)
3627 call getenv_loc('TORDPAR',tordname)
3628 open (itordp,file=tordname,status='old',readonly,shared)
3629 call getenv_loc('FOURIER',fouriername)
3630 open (ifourier,file=fouriername,status='old',readonly,shared)
3631 call getenv_loc('ELEPAR',elename)
3632 open (ielep,file=elename,status='old',readonly,shared)
3633 call getenv_loc('SIDEPAR',sidename)
3634 open (isidep,file=sidename,status='old',readonly,shared)
3635 #elif (defined CRAY) || (defined AIX)
3636 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3638 ! print *,"Processor",myrank," opened file 1"
3639 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3640 ! print *,"Processor",myrank," opened file 9"
3641 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3642 ! Get parameter filenames and open the parameter files.
3643 call getenv_loc('BONDPAR',bondname)
3644 open (ibond,file=bondname,status='old',action='read')
3645 ! print *,"Processor",myrank," opened file IBOND"
3646 call getenv_loc('THETPAR',thetname)
3647 open (ithep,file=thetname,status='old',action='read')
3648 ! print *,"Processor",myrank," opened file ITHEP"
3649 call getenv_loc('ROTPAR',rotname)
3650 open (irotam,file=rotname,status='old',action='read')
3651 ! print *,"Processor",myrank," opened file IROTAM"
3652 call getenv_loc('TORPAR',torname)
3653 open (itorp,file=torname,status='old',action='read')
3654 ! print *,"Processor",myrank," opened file ITORP"
3655 call getenv_loc('TORDPAR',tordname)
3656 open (itordp,file=tordname,status='old',action='read')
3657 ! print *,"Processor",myrank," opened file ITORDP"
3658 call getenv_loc('SCCORPAR',sccorname)
3659 open (isccor,file=sccorname,status='old',action='read')
3660 ! print *,"Processor",myrank," opened file ISCCOR"
3661 call getenv_loc('FOURIER',fouriername)
3662 open (ifourier,file=fouriername,status='old',action='read')
3663 ! print *,"Processor",myrank," opened file IFOURIER"
3664 call getenv_loc('ELEPAR',elename)
3665 open (ielep,file=elename,status='old',action='read')
3666 ! print *,"Processor",myrank," opened file IELEP"
3667 call getenv_loc('SIDEPAR',sidename)
3668 open (isidep,file=sidename,status='old',action='read')
3669 ! print *,"Processor",myrank," opened file ISIDEP"
3670 ! print *,"Processor",myrank," opened parameter files"
3672 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
3673 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3674 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3675 ! Get parameter filenames and open the parameter files.
3676 call getenv_loc('BONDPAR',bondname)
3677 open (ibond,file=bondname,status='old')
3678 call getenv_loc('THETPAR',thetname)
3679 open (ithep,file=thetname,status='old')
3680 call getenv_loc('ROTPAR',rotname)
3681 open (irotam,file=rotname,status='old')
3682 call getenv_loc('TORPAR',torname)
3683 open (itorp,file=torname,status='old')
3684 call getenv_loc('TORDPAR',tordname)
3685 open (itordp,file=tordname,status='old')
3686 call getenv_loc('SCCORPAR',sccorname)
3687 open (isccor,file=sccorname,status='old')
3688 call getenv_loc('FOURIER',fouriername)
3689 open (ifourier,file=fouriername,status='old')
3690 call getenv_loc('ELEPAR',elename)
3691 open (ielep,file=elename,status='old')
3692 call getenv_loc('SIDEPAR',sidename)
3693 open (isidep,file=sidename,status='old')
3695 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3697 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3698 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3699 ! Get parameter filenames and open the parameter files.
3700 call getenv_loc('BONDPAR',bondname)
3701 open (ibond,file=bondname,status='old',action='read')
3702 call getenv_loc('THETPAR',thetname)
3703 open (ithep,file=thetname,status='old',action='read')
3704 call getenv_loc('ROTPAR',rotname)
3705 open (irotam,file=rotname,status='old',action='read')
3706 call getenv_loc('TORPAR',torname)
3707 open (itorp,file=torname,status='old',action='read')
3708 call getenv_loc('TORDPAR',tordname)
3709 open (itordp,file=tordname,status='old',action='read')
3710 call getenv_loc('SCCORPAR',sccorname)
3711 open (isccor,file=sccorname,status='old',action='read')
3713 call getenv_loc('THETPARPDB',thetname_pdb)
3714 print *,"thetname_pdb ",thetname_pdb
3715 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
3716 print *,ithep_pdb," opened"
3718 call getenv_loc('FOURIER',fouriername)
3719 open (ifourier,file=fouriername,status='old',readonly)
3720 call getenv_loc('ELEPAR',elename)
3721 open (ielep,file=elename,status='old',readonly)
3722 call getenv_loc('SIDEPAR',sidename)
3723 open (isidep,file=sidename,status='old',readonly)
3725 call getenv_loc('ROTPARPDB',rotname_pdb)
3726 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
3731 ! 8/9/01 In the newest version SCp interaction constants are read from a file
3732 ! Use -DOLDSCP to use hard-coded constants instead.
3734 call getenv_loc('SCPPAR',scpname)
3735 #if defined(WINIFL) || defined(WINPGI)
3736 open (iscpp,file=scpname,status='old',readonly,shared)
3737 #elif (defined CRAY) || (defined AIX)
3738 open (iscpp,file=scpname,status='old',action='read')
3740 open (iscpp,file=scpname,status='old')
3742 open (iscpp,file=scpname,status='old',action='read')
3745 call getenv_loc('PATTERN',patname)
3746 #if defined(WINIFL) || defined(WINPGI)
3747 open (icbase,file=patname,status='old',readonly,shared)
3748 #elif (defined CRAY) || (defined AIX)
3749 open (icbase,file=patname,status='old',action='read')
3751 open (icbase,file=patname,status='old')
3753 open (icbase,file=patname,status='old',action='read')
3756 ! Open output file only for CG processes
3757 ! print *,"Processor",myrank," fg_rank",fg_rank
3758 if (fg_rank.eq.0) then
3760 if (nodes.eq.1) then
3763 npos = dlog10(dfloat(nodes-1))+1
3765 if (npos.lt.3) npos=3
3766 write (liczba,'(i1)') npos
3767 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
3769 write (liczba,form) me
3770 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
3771 liczba(:ilen(liczba))
3772 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3774 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3776 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
3777 liczba(:ilen(liczba))//'.mol2'
3778 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3779 liczba(:ilen(liczba))//'.stat'
3781 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
3782 //liczba(:ilen(liczba))//'.stat')
3783 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
3786 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3787 liczba(:ilen(liczba))//'.const'
3792 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
3793 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
3794 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
3795 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
3796 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
3798 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
3800 rest2name=prefix(:ilen(prefix))//'.rst'
3802 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
3805 #if defined(AIX) || defined(PGI)
3806 if (me.eq.king .or. .not. out1file) &
3807 open(iout,file=outname,status='unknown')
3809 if (fg_rank.gt.0) then
3810 write (liczba,'(i3.3)') myrank/nfgtasks
3811 write (ll,'(bz,i3.3)') fg_rank
3812 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3817 open(igeom,file=intname,status='unknown',position='append')
3818 open(ipdb,file=pdbname,status='unknown')
3819 open(imol2,file=mol2name,status='unknown')
3820 open(istat,file=statname,status='unknown',position='append')
3822 !1out open(iout,file=outname,status='unknown')
3825 if (me.eq.king .or. .not.out1file) &
3826 open(iout,file=outname,status='unknown')
3828 if (fg_rank.gt.0) then
3829 write (liczba,'(i3.3)') myrank/nfgtasks
3830 write (ll,'(bz,i3.3)') fg_rank
3831 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3836 open(igeom,file=intname,status='unknown',access='append')
3837 open(ipdb,file=pdbname,status='unknown')
3838 open(imol2,file=mol2name,status='unknown')
3839 open(istat,file=statname,status='unknown',access='append')
3841 !1out open(iout,file=outname,status='unknown')
3844 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
3845 csa_seed=prefix(:lenpre)//'.CSA.seed'
3846 csa_history=prefix(:lenpre)//'.CSA.history'
3847 csa_bank=prefix(:lenpre)//'.CSA.bank'
3848 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
3849 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
3850 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
3851 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
3852 csa_int=prefix(:lenpre)//'.int'
3853 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
3854 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
3855 csa_in=prefix(:lenpre)//'.CSA.in'
3856 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
3859 write (iout,'(80(1h-))')
3860 write (iout,'(30x,a)') "FILE ASSIGNMENT"
3861 write (iout,'(80(1h-))')
3862 write (iout,*) "Input file : ",&
3863 pref_orig(:ilen(pref_orig))//'.inp'
3864 write (iout,*) "Output file : ",&
3865 outname(:ilen(outname))
3867 write (iout,*) "Sidechain potential file : ",&
3868 sidename(:ilen(sidename))
3870 write (iout,*) "SCp potential file : ",&
3871 scpname(:ilen(scpname))
3873 write (iout,*) "Electrostatic potential file : ",&
3874 elename(:ilen(elename))
3875 write (iout,*) "Cumulant coefficient file : ",&
3876 fouriername(:ilen(fouriername))
3877 write (iout,*) "Torsional parameter file : ",&
3878 torname(:ilen(torname))
3879 write (iout,*) "Double torsional parameter file : ",&
3880 tordname(:ilen(tordname))
3881 write (iout,*) "SCCOR parameter file : ",&
3882 sccorname(:ilen(sccorname))
3883 write (iout,*) "Bond & inertia constant file : ",&
3884 bondname(:ilen(bondname))
3885 write (iout,*) "Bending parameter file : ",&
3886 thetname(:ilen(thetname))
3887 write (iout,*) "Rotamer parameter file : ",&
3888 rotname(:ilen(rotname))
3891 write (iout,*) "Thetpdb parameter file : ",&
3892 thetname_pdb(:ilen(thetname_pdb))
3895 write (iout,*) "Threading database : ",&
3896 patname(:ilen(patname))
3898 write (iout,*)" DIRTMP : ",&
3900 write (iout,'(80(1h-))')
3903 end subroutine openunits
3904 !-----------------------------------------------------------------------------
3907 use geometry_data, only: nres,dc
3909 ! implicit real*8 (a-h,o-z)
3910 ! include 'DIMENSIONS'
3911 ! include 'COMMON.CHAIN'
3912 ! include 'COMMON.IOUNITS'
3913 ! include 'COMMON.MD'
3916 ! real(kind=8) :: var,ene
3918 open(irest2,file=rest2name,status='unknown')
3919 read(irest2,*) totT,EK,potE,totE,t_bath
3921 ! AL 4/17/17: Now reading d_t(0,:) too
3923 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
3926 ! AL 4/17/17: Now reading d_c(0,:) too
3928 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
3931 read (irest2,*) iset
3935 end subroutine readrst
3936 !-----------------------------------------------------------------------------
3937 subroutine read_fragments
3941 use control_data, only:out1file
3944 ! implicit real*8 (a-h,o-z)
3945 ! include 'DIMENSIONS'
3949 ! include 'COMMON.SETUP'
3950 ! include 'COMMON.CHAIN'
3951 ! include 'COMMON.IOUNITS'
3952 ! include 'COMMON.MD'
3953 ! include 'COMMON.CONTROL'
3956 ! real(kind=8) :: var,ene
3958 read(inp,*) nset,nfrag,npair,nfrag_back
3960 !el from module energy
3961 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
3962 if(.not.allocated(wfrag_back)) then
3963 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
3964 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
3966 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
3967 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
3969 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
3970 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
3973 if(me.eq.king.or..not.out1file) &
3974 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
3975 " nfrag_back",nfrag_back
3977 read(inp,*) mset(iset)
3979 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
3981 if(me.eq.king.or..not.out1file) &
3982 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
3983 ifrag(2,i,iset), qinfrag(i,iset)
3986 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
3988 if(me.eq.king.or..not.out1file) &
3989 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
3990 ipair(2,i,iset), qinpair(i,iset)
3993 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
3994 wfrag_back(3,i,iset),&
3995 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
3996 if(me.eq.king.or..not.out1file) &
3997 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
3998 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4002 end subroutine read_fragments
4003 !-----------------------------------------------------------------------------
4005 !-----------------------------------------------------------------------------
4009 ! implicit real*8 (a-h,o-z)
4010 ! include 'DIMENSIONS'
4011 ! include 'COMMON.CSA'
4012 ! include 'COMMON.BANK'
4013 ! include 'COMMON.IOUNITS'
4015 ! integer :: ntf,ik,iw_pdb
4016 ! real(kind=8) :: var,ene
4018 open(icsa_in,file=csa_in,status="old",err=100)
4019 read(icsa_in,*) nconf
4020 read(icsa_in,*) jstart,jend
4021 read(icsa_in,*) nstmax
4022 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4023 read(icsa_in,*) nran0,nran1,irr
4024 read(icsa_in,*) nseed
4025 read(icsa_in,*) ntotal,cut1,cut2
4026 read(icsa_in,*) estop
4027 read(icsa_in,*) icmax,irestart
4028 read(icsa_in,*) ntbankm,dele,difcut
4029 read(icsa_in,*) iref,rmscut,pnccut
4030 read(icsa_in,*) ndiff
4037 end subroutine csa_read
4038 !-----------------------------------------------------------------------------
4039 subroutine initial_write
4042 ! implicit real*8 (a-h,o-z)
4043 ! include 'DIMENSIONS'
4044 ! include 'COMMON.CSA'
4045 ! include 'COMMON.BANK'
4046 ! include 'COMMON.IOUNITS'
4048 ! integer :: ntf,ik,iw_pdb
4049 ! real(kind=8) :: var,ene
4051 open(icsa_seed,file=csa_seed,status="unknown")
4052 write(icsa_seed,*) "seed"
4054 #if defined(AIX) || defined(PGI)
4055 open(icsa_history,file=csa_history,status="unknown",&
4058 open(icsa_history,file=csa_history,status="unknown",&
4061 write(icsa_history,*) nconf
4062 write(icsa_history,*) jstart,jend
4063 write(icsa_history,*) nstmax
4064 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4065 write(icsa_history,*) nran0,nran1,irr
4066 write(icsa_history,*) nseed
4067 write(icsa_history,*) ntotal,cut1,cut2
4068 write(icsa_history,*) estop
4069 write(icsa_history,*) icmax,irestart
4070 write(icsa_history,*) ntbankm,dele,difcut
4071 write(icsa_history,*) iref,rmscut,pnccut
4072 write(icsa_history,*) ndiff
4074 write(icsa_history,*)
4077 open(icsa_bank1,file=csa_bank1,status="unknown")
4078 write(icsa_bank1,*) 0
4082 end subroutine initial_write
4083 !-----------------------------------------------------------------------------
4084 subroutine restart_write
4087 ! implicit real*8 (a-h,o-z)
4088 ! include 'DIMENSIONS'
4089 ! include 'COMMON.IOUNITS'
4090 ! include 'COMMON.CSA'
4091 ! include 'COMMON.BANK'
4093 ! integer :: ntf,ik,iw_pdb
4094 ! real(kind=8) :: var,ene
4096 #if defined(AIX) || defined(PGI)
4097 open(icsa_history,file=csa_history,position="append")
4099 open(icsa_history,file=csa_history,access="append")
4101 write(icsa_history,*)
4102 write(icsa_history,*) "This is restart"
4103 write(icsa_history,*)
4104 write(icsa_history,*) nconf
4105 write(icsa_history,*) jstart,jend
4106 write(icsa_history,*) nstmax
4107 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4108 write(icsa_history,*) nran0,nran1,irr
4109 write(icsa_history,*) nseed
4110 write(icsa_history,*) ntotal,cut1,cut2
4111 write(icsa_history,*) estop
4112 write(icsa_history,*) icmax,irestart
4113 write(icsa_history,*) ntbankm,dele,difcut
4114 write(icsa_history,*) iref,rmscut,pnccut
4115 write(icsa_history,*) ndiff
4116 write(icsa_history,*)
4117 write(icsa_history,*) "irestart is: ", irestart
4119 write(icsa_history,*)
4123 end subroutine restart_write
4124 !-----------------------------------------------------------------------------
4126 !-----------------------------------------------------------------------------
4127 subroutine write_pdb(npdb,titelloc,ee)
4129 ! implicit real*8 (a-h,o-z)
4130 ! include 'DIMENSIONS'
4131 ! include 'COMMON.IOUNITS'
4132 character(len=50) :: titelloc1
4133 character*(*) :: titelloc
4134 character(len=3) :: zahl
4135 character(len=5) :: liczba5
4137 integer :: npdb !,ilen
4141 ! real(kind=8) :: var,ene
4145 if (npdb.lt.1000) then
4146 call numstr(npdb,zahl)
4147 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
4149 if (npdb.lt.10000) then
4150 write(liczba5,'(i1,i4)') 0,npdb
4152 write(liczba5,'(i5)') npdb
4154 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
4156 call pdbout(ee,titelloc1,ipdb)
4159 end subroutine write_pdb
4160 !-----------------------------------------------------------------------------
4162 !-----------------------------------------------------------------------------
4163 subroutine write_thread_summary
4164 ! Thread the sequence through a database of known structures
4165 use control_data, only: refstr
4167 use energy_data, only: n_ene_comp
4169 ! implicit real*8 (a-h,o-z)
4170 ! include 'DIMENSIONS'
4172 use MPI_data !include 'COMMON.INFO'
4175 ! include 'COMMON.CONTROL'
4176 ! include 'COMMON.CHAIN'
4177 ! include 'COMMON.DBASE'
4178 ! include 'COMMON.INTERACT'
4179 ! include 'COMMON.VAR'
4180 ! include 'COMMON.THREAD'
4181 ! include 'COMMON.FFIELD'
4182 ! include 'COMMON.SBRIDGE'
4183 ! include 'COMMON.HEADER'
4184 ! include 'COMMON.NAMES'
4185 ! include 'COMMON.IOUNITS'
4186 ! include 'COMMON.TIME1'
4188 integer,dimension(maxthread) :: ip
4189 real(kind=8),dimension(0:n_ene) :: energia
4191 integer :: i,j,ii,jj,ipj,ik,kk,ist
4192 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
4194 write (iout,'(30x,a/)') &
4195 ' *********** Summary threading statistics ************'
4196 write (iout,'(a)') 'Initial energies:'
4197 write (iout,'(a4,2x,a12,14a14,3a8)') &
4198 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
4199 'RMSnat','NatCONT','NNCONT','RMS'
4200 ! Energy sort patterns
4205 enet=ener(n_ene-1,ip(i))
4208 if (ener(n_ene-1,ip(j)).lt.enet) then
4210 enet=ener(n_ene-1,ip(j))
4222 ist=nres_base(2,ii)+ipatt(2,i)
4224 energia(i)=ener0(kk,i)
4226 etot=ener0(n_ene_comp+1,i)
4227 rmsnat=ener0(n_ene_comp+2,i)
4228 rms=ener0(n_ene_comp+3,i)
4229 frac=ener0(n_ene_comp+4,i)
4230 frac_nn=ener0(n_ene_comp+5,i)
4233 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4234 i,str_nam(ii),ist+1,&
4235 (energia(print_order(kk)),kk=1,nprint_ene),&
4236 etot,rmsnat,frac,frac_nn,rms
4238 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
4239 i,str_nam(ii),ist+1,&
4240 (energia(print_order(kk)),kk=1,nprint_ene),etot
4243 write (iout,'(//a)') 'Final energies:'
4244 write (iout,'(a4,2x,a12,17a14,3a8)') &
4245 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
4246 'RMSnat','NatCONT','NNCONT','RMS'
4250 ist=nres_base(2,ii)+ipatt(2,i)
4252 energia(kk)=ener(kk,ik)
4254 etot=ener(n_ene_comp+1,i)
4255 rmsnat=ener(n_ene_comp+2,i)
4256 rms=ener(n_ene_comp+3,i)
4257 frac=ener(n_ene_comp+4,i)
4258 frac_nn=ener(n_ene_comp+5,i)
4259 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4260 i,str_nam(ii),ist+1,&
4261 (energia(print_order(kk)),kk=1,nprint_ene),&
4262 etot,rmsnat,frac,frac_nn,rms
4264 write (iout,'(/a/)') 'IEXAM array:'
4265 write (iout,'(i5)') nexcl
4267 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
4269 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
4270 'Max. time for threading step ',max_time_for_thread,&
4271 'Average time for threading step: ',ave_time_for_thread
4273 end subroutine write_thread_summary
4274 !-----------------------------------------------------------------------------
4275 subroutine write_stat_thread(ithread,ipattern,ist)
4277 use energy_data, only: n_ene_comp
4279 ! implicit real*8 (a-h,o-z)
4280 ! include "DIMENSIONS"
4281 ! include "COMMON.CONTROL"
4282 ! include "COMMON.IOUNITS"
4283 ! include "COMMON.THREAD"
4284 ! include "COMMON.FFIELD"
4285 ! include "COMMON.DBASE"
4286 ! include "COMMON.NAMES"
4287 real(kind=8),dimension(0:n_ene) :: energia
4289 integer :: ithread,ipattern,ist,i
4290 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
4292 #if defined(AIX) || defined(PGI)
4293 open(istat,file=statname,position='append')
4295 open(istat,file=statname,access='append')
4298 energia(i)=ener(i,ithread)
4300 etot=ener(n_ene_comp+1,ithread)
4301 rmsnat=ener(n_ene_comp+2,ithread)
4302 rms=ener(n_ene_comp+3,ithread)
4303 frac=ener(n_ene_comp+4,ithread)
4304 frac_nn=ener(n_ene_comp+5,ithread)
4305 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4306 ithread,str_nam(ipattern),ist+1,&
4307 (energia(print_order(i)),i=1,nprint_ene),&
4308 etot,rmsnat,frac,frac_nn,rms
4311 end subroutine write_stat_thread
4312 !-----------------------------------------------------------------------------
4314 !-----------------------------------------------------------------------------
4315 end module io_config