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,1)),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,ncatprotparm
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,&
746 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
748 integer :: ichir1,ichir2
749 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
750 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
751 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
753 ! For printing parameters after they are read set the following in the UNRES
756 ! setenv PRINT_PARM YES
758 ! To print parameters in LaTeX format rather than as ASCII tables:
762 call getenv_loc("PRINT_PARM",lancuch)
763 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
764 call getenv_loc("LATEX",lancuch)
765 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 dwa16=2.0d0**(1.0d0/6.0d0)
769 ! Assign virtual-bond length
772 vblinv2=vblinv*vblinv
774 ! Read the virtual-bond parameters, masses, and moments of inertia
775 ! and Stokes' radii of the peptide group and side chains
777 allocate(dsc(ntyp1)) !(ntyp1)
778 allocate(dsc_inv(ntyp1)) !(ntyp1)
779 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
780 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
781 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
782 allocate(nbondterm(ntyp)) !(ntyp)
783 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
784 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
785 allocate(msc(ntyp+1,5)) !(ntyp+1)
786 allocate(isc(ntyp+1,5)) !(ntyp+1)
787 allocate(restok(ntyp+1,5)) !(ntyp+1)
788 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
789 allocate(long_r_sidechain(ntyp))
790 allocate(short_r_sidechain(ntyp))
795 read (ibond,*) vbldp0,akp,mp,ip,pstok
798 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
799 dsc(i) = vbldsc0(1,i)
803 dsc_inv(i)=1.0D0/dsc(i)
807 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
809 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
810 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
811 dsc(i) = vbldsc0(1,i)
815 dsc_inv(i)=1.0D0/dsc(i)
819 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
822 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
823 ! dsc(i) = vbldsc0_nucl(1,i)
827 ! dsc_inv(i)=1.0D0/dsc(i)
830 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
831 ! do i=1,ntyp_molec(2)
832 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
833 ! aksc_nucl(j,i),abond0_nucl(j,i),&
834 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
835 ! dsc(i) = vbldsc0(1,i)
839 ! dsc_inv(i)=1.0D0/dsc(i)
844 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
845 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
847 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
849 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
850 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
852 write (iout,'(13x,3f10.5)') &
853 vbldsc0(j,i),aksc(j,i),abond0(j,i)
858 read(iion,*) msc(i,5),restok(i,5)
859 print *,msc(i,5),restok(i,5)
863 read (iion,*) ncatprotparm
864 allocate(catprm(ncatprotparm,4))
866 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
869 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
870 !----------------------------------------------------
871 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
872 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
873 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
874 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
875 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
876 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
886 allocate(liptranene(ntyp))
887 !C reading lipid parameters
888 write (iout,*) "iliptranpar",iliptranpar
890 read(iliptranpar,*) pepliptran
893 read(iliptranpar,*) liptranene(i)
894 print *,liptranene(i)
900 ! Read the parameters of the probability distribution/energy expression
901 ! of the virtual-bond valence angles theta
904 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
905 (bthet(j,i,1,1),j=1,2)
906 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
907 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
908 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
912 athet(1,i,1,-1)=athet(1,i,1,1)
913 athet(2,i,1,-1)=athet(2,i,1,1)
914 bthet(1,i,1,-1)=-bthet(1,i,1,1)
915 bthet(2,i,1,-1)=-bthet(2,i,1,1)
916 athet(1,i,-1,1)=-athet(1,i,1,1)
917 athet(2,i,-1,1)=-athet(2,i,1,1)
918 bthet(1,i,-1,1)=bthet(1,i,1,1)
919 bthet(2,i,-1,1)=bthet(2,i,1,1)
923 athet(1,i,-1,-1)=athet(1,-i,1,1)
924 athet(2,i,-1,-1)=-athet(2,-i,1,1)
925 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
926 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
927 athet(1,i,-1,1)=athet(1,-i,1,1)
928 athet(2,i,-1,1)=-athet(2,-i,1,1)
929 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
930 bthet(2,i,-1,1)=bthet(2,-i,1,1)
931 athet(1,i,1,-1)=-athet(1,-i,1,1)
932 athet(2,i,1,-1)=athet(2,-i,1,1)
933 bthet(1,i,1,-1)=bthet(1,-i,1,1)
934 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
939 polthet(j,i)=polthet(j,-i)
942 gthet(j,i)=gthet(j,-i)
950 'Parameters of the virtual-bond valence angles:'
951 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
952 ' ATHETA0 ',' A1 ',' A2 ',&
955 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
956 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
958 write (iout,'(/a/9x,5a/79(1h-))') &
959 'Parameters of the expression for sigma(theta_c):',&
960 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
961 ' ALPH3 ',' SIGMA0C '
963 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
964 (polthet(j,i),j=0,3),sigc0(i)
966 write (iout,'(/a/9x,5a/79(1h-))') &
967 'Parameters of the second gaussian:',&
968 ' THETA0 ',' SIGMA0 ',' G1 ',&
971 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
972 sig0(i),(gthet(j,i),j=1,3)
976 'Parameters of the virtual-bond valence angles:'
977 write (iout,'(/a/9x,5a/79(1h-))') &
978 'Coefficients of expansion',&
979 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
980 ' b1*10^1 ',' b2*10^1 '
982 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
983 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
984 (10*bthet(j,i,1,1),j=1,2)
986 write (iout,'(/a/9x,5a/79(1h-))') &
987 'Parameters of the expression for sigma(theta_c):',&
988 ' alpha0 ',' alph1 ',' alph2 ',&
989 ' alhp3 ',' sigma0c '
991 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
992 (polthet(j,i),j=0,3),sigc0(i)
994 write (iout,'(/a/9x,5a/79(1h-))') &
995 'Parameters of the second gaussian:',&
996 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
999 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1000 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1006 ! Read the parameters of Utheta determined from ab initio surfaces
1007 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1009 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1010 ntheterm3,nsingle,ndouble
1011 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1013 !----------------------------------------------------
1014 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1015 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1016 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1017 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1018 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1019 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1020 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1021 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1022 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1023 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1024 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1025 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1026 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1027 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1028 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1029 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1030 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1031 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1032 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1033 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1034 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1035 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1036 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1037 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1039 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1041 ithetyp(i)=-ithetyp(-i)
1044 aa0thet(:,:,:,:)=0.0d0
1045 aathet(:,:,:,:,:)=0.0d0
1046 bbthet(:,:,:,:,:,:)=0.0d0
1047 ccthet(:,:,:,:,:,:)=0.0d0
1048 ddthet(:,:,:,:,:,:)=0.0d0
1049 eethet(:,:,:,:,:,:)=0.0d0
1050 ffthet(:,:,:,:,:,:,:)=0.0d0
1051 ggthet(:,:,:,:,:,:,:)=0.0d0
1053 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1055 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1056 ! VAR:1=non-glicyne non-proline 2=proline
1057 ! VAR:negative values for D-aminoacid
1059 do j=-nthetyp,nthetyp
1060 do k=-nthetyp,nthetyp
1061 read (ithep,'(6a)',end=111,err=111) res1
1062 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1063 ! VAR: aa0thet is variable describing the average value of Foureir
1064 ! VAR: expansion series
1065 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1066 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1067 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1068 read (ithep,*,end=111,err=111) &
1069 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1070 read (ithep,*,end=111,err=111) &
1071 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1072 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1073 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1074 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1076 read (ithep,*,end=111,err=111) &
1077 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1078 ffthet(lll,llll,ll,i,j,k,iblock),&
1079 ggthet(llll,lll,ll,i,j,k,iblock),&
1080 ggthet(lll,llll,ll,i,j,k,iblock),&
1081 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1086 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1087 ! coefficients of theta-and-gamma-dependent terms are zero.
1088 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1089 ! RECOMENTDED AFTER VERSION 3.3)
1093 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1094 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1096 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1097 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1100 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1102 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1105 ! AND COMMENT THE LOOPS BELOW
1109 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1110 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1112 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1113 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1116 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1118 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1123 ! Substitution for D aminoacids from symmetry.
1126 do j=-nthetyp,nthetyp
1127 do k=-nthetyp,nthetyp
1128 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1130 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1134 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1135 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1136 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1137 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1143 ffthet(llll,lll,ll,i,j,k,iblock)= &
1144 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1145 ffthet(lll,llll,ll,i,j,k,iblock)= &
1146 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1147 ggthet(llll,lll,ll,i,j,k,iblock)= &
1148 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1149 ggthet(lll,llll,ll,i,j,k,iblock)= &
1150 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1159 ! Control printout of the coefficients of virtual-bond-angle potentials
1162 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1167 write (iout,'(//4a)') &
1168 'Type ',onelett(i),onelett(j),onelett(k)
1169 write (iout,'(//a,10x,a)') " l","a[l]"
1170 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1171 write (iout,'(i2,1pe15.5)') &
1172 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1174 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1175 "b",l,"c",l,"d",l,"e",l
1177 write (iout,'(i2,4(1pe15.5))') m,&
1178 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1179 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1183 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1184 "f+",l,"f-",l,"g+",l,"g-",l
1187 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1188 ffthet(n,m,l,i,j,k,iblock),&
1189 ffthet(m,n,l,i,j,k,iblock),&
1190 ggthet(n,m,l,i,j,k,iblock),&
1191 ggthet(m,n,l,i,j,k,iblock)
1201 write (2,*) "Start reading THETA_PDB",ithep_pdb
1203 ! write (2,*) 'i=',i
1204 read (ithep_pdb,*,err=111,end=111) &
1205 a0thet(i),(athet(j,i,1,1),j=1,2),&
1206 (bthet(j,i,1,1),j=1,2)
1207 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1208 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1209 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1210 sigc0(i)=sigc0(i)**2
1213 athet(1,i,1,-1)=athet(1,i,1,1)
1214 athet(2,i,1,-1)=athet(2,i,1,1)
1215 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1216 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1217 athet(1,i,-1,1)=-athet(1,i,1,1)
1218 athet(2,i,-1,1)=-athet(2,i,1,1)
1219 bthet(1,i,-1,1)=bthet(1,i,1,1)
1220 bthet(2,i,-1,1)=bthet(2,i,1,1)
1223 a0thet(i)=a0thet(-i)
1224 athet(1,i,-1,-1)=athet(1,-i,1,1)
1225 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1226 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1227 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1228 athet(1,i,-1,1)=athet(1,-i,1,1)
1229 athet(2,i,-1,1)=-athet(2,-i,1,1)
1230 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1231 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1232 athet(1,i,1,-1)=-athet(1,-i,1,1)
1233 athet(2,i,1,-1)=athet(2,-i,1,1)
1234 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1235 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1236 theta0(i)=theta0(-i)
1240 polthet(j,i)=polthet(j,-i)
1243 gthet(j,i)=gthet(j,-i)
1246 write (2,*) "End reading THETA_PDB"
1250 !--------------- Reading theta parameters for nucleic acid-------
1251 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1252 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1253 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1254 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1255 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1256 nthetyp_nucl+1,nthetyp_nucl+1))
1257 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1258 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1259 nthetyp_nucl+1,nthetyp_nucl+1))
1260 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1261 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1262 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1263 nthetyp_nucl+1,nthetyp_nucl+1))
1264 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1265 nthetyp_nucl+1,nthetyp_nucl+1))
1266 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1267 nthetyp_nucl+1,nthetyp_nucl+1))
1268 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1269 nthetyp_nucl+1,nthetyp_nucl+1))
1270 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1271 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1272 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1273 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1274 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1275 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1277 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1278 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1280 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1282 aa0thet_nucl(:,:,:)=0.0d0
1283 aathet_nucl(:,:,:,:)=0.0d0
1284 bbthet_nucl(:,:,:,:,:)=0.0d0
1285 ccthet_nucl(:,:,:,:,:)=0.0d0
1286 ddthet_nucl(:,:,:,:,:)=0.0d0
1287 eethet_nucl(:,:,:,:,:)=0.0d0
1288 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1289 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1294 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1295 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1296 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1297 read (ithep_nucl,*,end=111,err=111) &
1298 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1299 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1300 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1301 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1302 read (ithep_nucl,*,end=111,err=111) &
1303 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1304 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1305 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1310 !-------------------------------------------
1311 allocate(nlob(ntyp1)) !(ntyp1)
1312 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1313 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1314 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1321 gaussc(:,:,:,:)=0.0D0
1325 ! Read the parameters of the probability distribution/energy expression
1326 ! of the side chains.
1329 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1333 dsc_inv(i)=1.0D0/dsc(i)
1344 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1345 ((blower(k,l,1),l=1,k),k=1,3)
1346 censc(1,1,-i)=censc(1,1,i)
1347 censc(2,1,-i)=censc(2,1,i)
1348 censc(3,1,-i)=-censc(3,1,i)
1350 read (irotam,*,end=112,err=112) bsc(j,i)
1351 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1352 ((blower(k,l,j),l=1,k),k=1,3)
1353 censc(1,j,-i)=censc(1,j,i)
1354 censc(2,j,-i)=censc(2,j,i)
1355 censc(3,j,-i)=-censc(3,j,i)
1356 ! BSC is amplitude of Gaussian
1363 akl=akl+blower(k,m,j)*blower(l,m,j)
1367 if (((k.eq.3).and.(l.ne.3)) &
1368 .or.((l.eq.3).and.(k.ne.3))) then
1369 gaussc(k,l,j,-i)=-akl
1370 gaussc(l,k,j,-i)=-akl
1372 gaussc(k,l,j,-i)=akl
1373 gaussc(l,k,j,-i)=akl
1382 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1385 if (nlobi.gt.0) then
1387 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1388 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1389 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1390 'log h',(bsc(j,i),j=1,nlobi)
1391 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1392 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1394 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1395 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1398 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1399 write (iout,'(a,f10.4,4(16x,f10.4))') &
1400 'Center ',(bsc(j,i),j=1,nlobi)
1401 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1410 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1411 ! added by Urszula Kozlowska 07/11/2007
1413 !el Maximum number of SC local term fitting function coefficiants
1414 !el integer,parameter :: maxsccoef=65
1416 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1419 read (irotam,*,end=112,err=112)
1421 read (irotam,*,end=112,err=112)
1424 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1428 !---------reading nucleic acid parameters for rotamers-------------------
1429 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1430 do i=1,ntyp_molec(2)
1431 read (irotam_nucl,*,end=112,err=112)
1433 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1439 write (iout,*) "Base rotamer parameters"
1440 do i=1,ntyp_molec(2)
1441 write (iout,'(a)') restyp(i,2)
1442 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1447 ! Read the parameters of the probability distribution/energy expression
1448 ! of the side chains.
1450 write (2,*) "Start reading ROTAM_PDB"
1452 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1456 dsc_inv(i)=1.0D0/dsc(i)
1467 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1468 ((blower(k,l,1),l=1,k),k=1,3)
1470 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1471 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1472 ((blower(k,l,j),l=1,k),k=1,3)
1479 akl=akl+blower(k,m,j)*blower(l,m,j)
1489 write (2,*) "End reading ROTAM_PDB"
1495 ! Read torsional parameters in old format
1497 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1499 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1500 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1501 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1503 !el from energy module--------
1504 allocate(v1(nterm_old,ntortyp,ntortyp))
1505 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1506 !el---------------------------
1511 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1517 write (iout,'(/a/)') 'Torsional constants:'
1520 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1521 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1527 ! Read torsional parameters
1529 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1530 read (itorp,*,end=113,err=113) ntortyp
1531 !el from energy module---------
1532 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1533 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1535 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1536 allocate(vlor2(maxlor,ntortyp,ntortyp))
1537 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1538 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1540 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1541 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1542 !el---------------------------
1545 !el---------------------------
1547 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1549 itortyp(i)=-itortyp(-i)
1551 itortyp(ntyp1)=ntortyp
1552 itortyp(-ntyp1)=-ntortyp
1554 write (iout,*) 'ntortyp',ntortyp
1556 do j=-ntortyp+1,ntortyp-1
1557 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1559 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1560 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1563 do k=1,nterm(i,j,iblock)
1564 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1566 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1567 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1568 v0ij=v0ij+si*v1(k,i,j,iblock)
1570 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1571 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1572 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1574 do k=1,nlor(i,j,iblock)
1575 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1576 vlor2(k,i,j),vlor3(k,i,j)
1577 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1580 v0(-i,-j,iblock)=v0ij
1586 write (iout,'(/a/)') 'Torsional constants:'
1588 do i=-ntortyp,ntortyp
1589 do j=-ntortyp,ntortyp
1590 write (iout,*) 'ityp',i,' jtyp',j
1591 write (iout,*) 'Fourier constants'
1592 do k=1,nterm(i,j,iblock)
1593 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1596 write (iout,*) 'Lorenz constants'
1597 do k=1,nlor(i,j,iblock)
1598 write (iout,'(3(1pe15.5))') &
1599 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1605 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1607 ! 6/23/01 Read parameters for double torsionals
1609 !el from energy module------------
1610 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1611 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1612 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1613 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1614 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1615 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1616 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1617 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1618 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1619 !---------------------------------
1623 do j=-ntortyp+1,ntortyp-1
1624 do k=-ntortyp+1,ntortyp-1
1625 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1626 ! write (iout,*) "OK onelett",
1629 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1630 .or. t3.ne.toronelet(k)) then
1631 write (iout,*) "Error in double torsional parameter file",&
1634 call MPI_Finalize(Ierror)
1636 stop "Error in double torsional parameter file"
1638 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1639 ntermd_2(i,j,k,iblock)
1640 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1641 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1642 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1643 ntermd_1(i,j,k,iblock))
1644 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1645 ntermd_1(i,j,k,iblock))
1646 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1647 ntermd_1(i,j,k,iblock))
1648 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1649 ntermd_1(i,j,k,iblock))
1650 ! Martix of D parameters for one dimesional foureir series
1651 do l=1,ntermd_1(i,j,k,iblock)
1652 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1653 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1654 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1655 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1656 ! write(iout,*) "whcodze" ,
1657 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1659 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1660 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1661 v2s(m,l,i,j,k,iblock),&
1662 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1663 ! Martix of D parameters for two dimesional fourier series
1664 do l=1,ntermd_2(i,j,k,iblock)
1666 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1667 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1668 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1669 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1678 write (iout,*) 'Constants for double torsionals'
1681 do j=-ntortyp+1,ntortyp-1
1682 do k=-ntortyp+1,ntortyp-1
1683 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1684 ' nsingle',ntermd_1(i,j,k,iblock),&
1685 ' ndouble',ntermd_2(i,j,k,iblock)
1687 write (iout,*) 'Single angles:'
1688 do l=1,ntermd_1(i,j,k,iblock)
1689 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1690 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1691 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1692 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1695 write (iout,*) 'Pairs of angles:'
1696 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1697 do l=1,ntermd_2(i,j,k,iblock)
1698 write (iout,'(i5,20f10.5)') &
1699 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1702 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1703 do l=1,ntermd_2(i,j,k,iblock)
1704 write (iout,'(i5,20f10.5)') &
1705 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1706 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1715 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1716 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
1717 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
1718 !el from energy module---------
1719 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1720 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1722 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1723 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
1724 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
1725 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1727 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
1728 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1729 !el---------------------------
1732 !el--------------------
1733 read (itorp_nucl,*,end=113,err=113) &
1734 (itortyp_nucl(i),i=1,ntyp_molec(2))
1735 ! print *,itortyp_nucl(:)
1736 !c write (iout,*) 'ntortyp',ntortyp
1739 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
1740 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
1743 do k=1,nterm_nucl(i,j)
1744 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
1745 v0ij=v0ij+si*v1_nucl(k,i,j)
1748 do k=1,nlor_nucl(i,j)
1749 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
1750 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
1751 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
1757 ! Read of Side-chain backbone correlation parameters
1758 ! Modified 11 May 2012 by Adasko
1761 read (isccor,*,end=119,err=119) nsccortyp
1763 !el from module energy-------------
1764 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1765 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1766 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1767 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1768 !-----------------------------------
1770 !el from module energy-------------
1771 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1773 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1775 isccortyp(i)=-isccortyp(-i)
1777 iscprol=isccortyp(20)
1778 ! write (iout,*) 'ntortyp',ntortyp
1780 !c maxinter is maximum interaction sites
1781 !el from module energy---------
1782 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1783 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1784 -nsccortyp:nsccortyp))
1785 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1786 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1787 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1788 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1789 !-----------------------------------
1793 read (isccor,*,end=119,err=119) &
1794 nterm_sccor(i,j),nlor_sccor(i,j)
1800 nterm_sccor(-i,j)=nterm_sccor(i,j)
1801 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1802 nterm_sccor(i,-j)=nterm_sccor(i,j)
1803 do k=1,nterm_sccor(i,j)
1804 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1806 if (j.eq.iscprol) then
1807 if (i.eq.isccortyp(10)) then
1808 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1809 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1811 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1812 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1813 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1814 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1815 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1816 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1817 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1818 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1821 if (i.eq.isccortyp(10)) then
1822 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1823 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1825 if (j.eq.isccortyp(10)) then
1826 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1827 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1829 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1830 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1831 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1832 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1833 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1834 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1838 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1839 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1840 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1841 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1844 do k=1,nlor_sccor(i,j)
1845 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1846 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1847 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1848 (1+vlor3sccor(k,i,j)**2)
1850 v0sccor(l,i,j)=v0ijsccor
1851 v0sccor(l,-i,j)=v0ijsccor1
1852 v0sccor(l,i,-j)=v0ijsccor2
1853 v0sccor(l,-i,-j)=v0ijsccor3
1859 !el from module energy-------------
1860 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1862 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1863 ! write (iout,*) 'ntortyp',ntortyp
1865 !c maxinter is maximum interaction sites
1866 !el from module energy---------
1867 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1868 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1869 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1870 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1871 !-----------------------------------
1875 read (isccor,*,end=119,err=119) &
1876 nterm_sccor(i,j),nlor_sccor(i,j)
1880 do k=1,nterm_sccor(i,j)
1881 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1883 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1886 do k=1,nlor_sccor(i,j)
1887 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1888 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1889 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1890 (1+vlor3sccor(k,i,j)**2)
1892 v0sccor(l,i,j)=v0ijsccor !el ,iblock
1900 write (iout,'(/a/)') 'Torsional constants:'
1903 write (iout,*) 'ityp',i,' jtyp',j
1904 write (iout,*) 'Fourier constants'
1905 do k=1,nterm_sccor(i,j)
1906 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1908 write (iout,*) 'Lorenz constants'
1909 do k=1,nlor_sccor(i,j)
1910 write (iout,'(3(1pe15.5))') &
1911 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1918 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1919 ! interaction energy of the Gly, Ala, and Pro prototypes.
1924 write (iout,*) "Coefficients of the cumulants"
1926 read (ifourier,*) nloctyp
1927 !write(iout,*) "nloctyp",nloctyp
1928 !el from module energy-------
1929 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1930 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1931 allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1932 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1933 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1934 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1935 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1936 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1940 !--------------------------------
1943 read (ifourier,*,end=115,err=115)
1944 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
1946 write (iout,*) 'Type',i
1947 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
1957 B1tilde(1,-i) =-b(3)
1959 ! b1tilde(1,i)=0.0d0
1960 ! b1tilde(2,i)=0.0d0
1985 Ctilde(1,2,-i)=-b(9)
1989 ! Ctilde(1,1,i)=0.0d0
1990 ! Ctilde(1,2,i)=0.0d0
1991 ! Ctilde(2,1,i)=0.0d0
1992 ! Ctilde(2,2,i)=0.0d0
2010 Dtilde(1,2,-i)=-b(8)
2014 ! Dtilde(1,1,i)=0.0d0
2015 ! Dtilde(1,2,i)=0.0d0
2016 ! Dtilde(2,1,i)=0.0d0
2017 ! Dtilde(2,2,i)=0.0d0
2018 EE(1,1,i)= b(10)+b(11)
2019 EE(2,2,i)=-b(10)+b(11)
2020 EE(2,1,i)= b(12)-b(13)
2021 EE(1,2,i)= b(12)+b(13)
2022 EE(1,1,-i)= b(10)+b(11)
2023 EE(2,2,-i)=-b(10)+b(11)
2024 EE(2,1,-i)=-b(12)+b(13)
2025 EE(1,2,-i)=-b(12)-b(13)
2031 ! ee(2,1,i)=ee(1,2,i)
2035 write (iout,*) 'Type',i
2037 write(iout,*) B1(1,i),B1(2,i)
2039 write(iout,*) B2(1,i),B2(2,i)
2042 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
2046 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
2050 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
2055 ! Read electrostatic-interaction parameters
2060 write (iout,'(/a)') 'Electrostatic interaction constants:'
2061 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2062 'IT','JT','APP','BPP','AEL6','AEL3'
2064 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2065 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2066 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2067 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2072 app (i,j)=epp(i,j)*rri*rri
2073 bpp (i,j)=-2.0D0*epp(i,j)*rri
2074 ael6(i,j)=elpp6(i,j)*4.2D0**6
2075 ael3(i,j)=elpp3(i,j)*4.2D0**3
2077 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2083 ! Read side-chain interaction parameters.
2085 !el from module energy - COMMON.INTERACT-------
2086 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2087 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2088 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2090 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2091 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2092 allocate(epslip(ntyp,ntyp))
2100 !--------------------------------
2102 read (isidep,*,end=117,err=117) ipot,expon
2103 if (ipot.lt.1 .or. ipot.gt.5) then
2104 write (iout,'(2a)') 'Error while reading SC interaction',&
2105 'potential file - unknown potential type.'
2107 call MPI_Finalize(Ierror)
2113 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2114 ', exponents are ',expon,2*expon
2115 ! goto (10,20,30,30,40) ipot
2117 !----------------------- LJ potential ---------------------------------
2119 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2120 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2121 (sigma0(i),i=1,ntyp)
2123 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2124 write (iout,'(a/)') 'The epsilon array:'
2125 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2126 write (iout,'(/a)') 'One-body parameters:'
2127 write (iout,'(a,4x,a)') 'residue','sigma'
2128 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2131 !----------------------- LJK potential --------------------------------
2133 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2134 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2135 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2137 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2138 write (iout,'(a/)') 'The epsilon array:'
2139 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2140 write (iout,'(/a)') 'One-body parameters:'
2141 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2142 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2146 !---------------------- GB or BP potential -----------------------------
2150 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2152 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2153 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2154 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2155 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2157 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2160 ! For the GB potential convert sigma'**2 into chi'
2163 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2167 write (iout,'(/a/)') 'Parameters of the BP potential:'
2168 write (iout,'(a/)') 'The epsilon array:'
2169 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2170 write (iout,'(/a)') 'One-body parameters:'
2171 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2173 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2174 chip(i),alp(i),i=1,ntyp)
2177 !--------------------- GBV potential -----------------------------------
2179 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2180 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2181 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2182 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2184 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2185 write (iout,'(a/)') 'The epsilon array:'
2186 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2187 write (iout,'(/a)') 'One-body parameters:'
2188 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2189 's||/s_|_^2',' chip ',' alph '
2190 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2191 sigii(i),chip(i),alp(i),i=1,ntyp)
2194 write(iout,*)"Wrong ipot"
2200 !-----------------------------------------------------------------------
2201 ! Calculate the "working" parameters of SC interactions.
2203 !el from module energy - COMMON.INTERACT-------
2204 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2205 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2206 allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2207 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2209 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2222 sc_aa_tube_par(:)=0.0d0
2223 sc_bb_tube_par(:)=0.0d0
2225 !--------------------------------
2230 epslip(i,j)=epslip(j,i)
2235 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2236 sigma(j,i)=sigma(i,j)
2237 rs0(i,j)=dwa16*sigma(i,j)
2241 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2242 'Working parameters of the SC interactions:',&
2243 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2248 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2257 sigeps=dsign(1.0D0,epsij)
2259 aa_aq(i,j)=epsij*rrij*rrij
2260 bb_aq(i,j)=-sigeps*epsij*rrij
2261 aa_aq(j,i)=aa_aq(i,j)
2262 bb_aq(j,i)=bb_aq(i,j)
2263 epsijlip=epslip(i,j)
2264 sigeps=dsign(1.0D0,epsijlip)
2265 epsijlip=dabs(epsijlip)
2266 aa_lip(i,j)=epsijlip*rrij*rrij
2267 bb_lip(i,j)=-sigeps*epsijlip*rrij
2268 aa_lip(j,i)=aa_lip(i,j)
2269 bb_lip(j,i)=bb_lip(i,j)
2270 !C write(iout,*) aa_lip
2272 sigt1sq=sigma0(i)**2
2273 sigt2sq=sigma0(j)**2
2276 ratsig1=sigt2sq/sigt1sq
2277 ratsig2=1.0D0/ratsig1
2278 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2279 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2280 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2284 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2285 sigmaii(i,j)=rsum_max
2286 sigmaii(j,i)=rsum_max
2288 ! sigmaii(i,j)=r0(i,j)
2289 ! sigmaii(j,i)=r0(i,j)
2291 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2292 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2293 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2294 augm(i,j)=epsij*r_augm**(2*expon)
2295 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2302 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2303 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2304 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2309 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2310 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2311 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2312 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2313 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2314 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2315 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2316 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2317 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2318 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2319 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2320 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2321 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2322 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2323 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2324 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2333 read (isidep_nucl,*) ipot_nucl
2334 ! print *,"TU?!",ipot_nucl
2335 if (ipot_nucl.eq.1) then
2336 do i=1,ntyp_molec(2)
2337 do j=i,ntyp_molec(2)
2338 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2339 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2343 do i=1,ntyp_molec(2)
2344 do j=i,ntyp_molec(2)
2345 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2346 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2347 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2351 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2352 do i=1,ntyp_molec(2)
2353 do j=i,ntyp_molec(2)
2354 rrij=sigma_nucl(i,j)
2358 epsij=4*eps_nucl(i,j)
2359 sigeps=dsign(1.0D0,epsij)
2361 aa_nucl(i,j)=epsij*rrij*rrij
2362 bb_nucl(i,j)=-sigeps*epsij*rrij
2363 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2364 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2365 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2366 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2367 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2368 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2371 aa_nucl(i,j)=aa_nucl(j,i)
2372 bb_nucl(i,j)=bb_nucl(j,i)
2373 ael3_nucl(i,j)=ael3_nucl(j,i)
2374 ael6_nucl(i,j)=ael6_nucl(j,i)
2375 ael63_nucl(i,j)=ael63_nucl(j,i)
2376 ael32_nucl(i,j)=ael32_nucl(j,i)
2377 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2378 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2379 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2380 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2381 eps_nucl(i,j)=eps_nucl(j,i)
2382 sigma_nucl(i,j)=sigma_nucl(j,i)
2383 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2387 write(iout,*) "tube param"
2388 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2389 ccavtubpep,dcavtubpep,tubetranenepep
2390 sigmapeptube=sigmapeptube**6
2391 sigeps=dsign(1.0D0,epspeptube)
2392 epspeptube=dabs(epspeptube)
2393 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2394 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2395 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2397 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2398 ccavtub(i),dcavtub(i),tubetranene(i)
2399 sigmasctube=sigmasctube**6
2400 sigeps=dsign(1.0D0,epssctube)
2401 epssctube=dabs(epssctube)
2402 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2403 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2404 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2406 !-----------------READING SC BASE POTENTIALS-----------------------------
2407 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2408 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2409 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2410 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2411 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2412 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2413 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2414 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2415 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2416 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2417 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2418 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2419 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2420 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2421 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2422 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2425 do i=1,ntyp_molec(1)
2426 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2427 write (*,*) "Im in ", i, " ", j
2428 read(isidep_scbase,*) &
2429 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2430 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2431 write(*,*) "eps",eps_scbase(i,j)
2432 read(isidep_scbase,*) &
2433 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2434 chis_scbase(i,j,1),chis_scbase(i,j,2)
2435 read(isidep_scbase,*) &
2436 dhead_scbasei(i,j), &
2437 dhead_scbasej(i,j), &
2438 rborn_scbasei(i,j),rborn_scbasej(j,i)
2439 read(isidep_scbase,*) &
2440 (wdipdip_scbase(k,i,j),k=1,3), &
2441 (wqdip_scbase(k,i,j),k=1,2)
2442 read(isidep_scbase,*) &
2443 alphapol_scbase(i,j), &
2444 epsintab_scbase(i,j)
2447 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2448 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2450 do i=1,ntyp_molec(1)
2451 do j=1,ntyp_molec(2)-1
2452 epsij=eps_scbase(i,j)
2458 sigeps=dsign(1.0D0,epsij)
2460 aa_scbase(i,j)=epsij*rrij*rrij
2461 bb_scbase(i,j)=-sigeps*epsij*rrij
2465 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2470 ! Define the SC-p interaction constants (hard-coded; old style)
2473 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2475 ! aad(i,1)=0.3D0*4.0D0**12
2476 ! Following line for constants currently implemented
2477 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2478 aad(i,1)=1.5D0*4.0D0**12
2479 ! aad(i,1)=0.17D0*5.6D0**12
2481 ! "Soft" SC-p repulsion
2483 ! Following line for constants currently implemented
2484 ! aad(i,1)=0.3D0*4.0D0**6
2485 ! "Hard" SC-p repulsion
2486 bad(i,1)=3.0D0*4.0D0**6
2487 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2496 ! 8/9/01 Read the SC-p interaction constants from file
2499 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2502 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2503 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2504 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2505 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2509 write (iout,*) "Parameters of SC-p interactions:"
2511 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2512 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2517 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
2519 do i=1,ntyp_molec(2)
2520 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
2522 do i=1,ntyp_molec(2)
2523 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
2524 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
2526 r0pp=1.12246204830937298142*5.16158
2532 ! Define the constants of the disulfide bridge
2536 ! Old arbitrary potential - commented out.
2541 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2542 ! energy surface of diethyl disulfide.
2543 ! A. Liwo and U. Kozlowska, 11/24/03
2560 write (iout,'(/a)') "Disulfide bridge parameters:"
2561 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2562 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2563 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2564 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2568 111 write (iout,*) "Error reading bending energy parameters."
2570 112 write (iout,*) "Error reading rotamer energy parameters."
2572 113 write (iout,*) "Error reading torsional energy parameters."
2574 114 write (iout,*) "Error reading double torsional energy parameters."
2576 115 write (iout,*) &
2577 "Error reading cumulant (multibody energy) parameters."
2579 116 write (iout,*) "Error reading electrostatic energy parameters."
2581 117 write (iout,*) "Error reading side chain interaction parameters."
2583 118 write (iout,*) "Error reading SCp interaction parameters."
2585 119 write (iout,*) "Error reading SCCOR parameters"
2588 call MPI_Finalize(Ierror)
2592 end subroutine parmread
2594 !-----------------------------------------------------------------------------
2596 !-----------------------------------------------------------------------------
2597 subroutine printmat(ldim,m,n,iout,key,a)
2600 character(len=3),dimension(n) :: key
2601 real(kind=8),dimension(ldim,n) :: a
2603 integer :: i,j,k,m,iout,nlim
2607 write (iout,1000) (key(k),k=i,nlim)
2609 1000 format (/5x,8(6x,a3))
2610 1020 format (/80(1h-)/)
2612 write (iout,1010) key(j),(a(j,k),k=i,nlim)
2615 1010 format (a3,2x,8(f9.4))
2617 end subroutine printmat
2618 !-----------------------------------------------------------------------------
2620 !-----------------------------------------------------------------------------
2622 ! Read the PDB file and convert the peptide geometry into virtual-chain
2625 use energy_data, only: itype,istype
2629 use control, only: rescode,sugarcode
2630 ! implicit real*8 (a-h,o-z)
2631 ! include 'DIMENSIONS'
2632 ! include 'COMMON.LOCAL'
2633 ! include 'COMMON.VAR'
2634 ! include 'COMMON.CHAIN'
2635 ! include 'COMMON.INTERACT'
2636 ! include 'COMMON.IOUNITS'
2637 ! include 'COMMON.GEO'
2638 ! include 'COMMON.NAMES'
2639 ! include 'COMMON.CONTROL'
2640 ! include 'COMMON.DISTFIT'
2641 ! include 'COMMON.SETUP'
2642 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
2644 logical :: lprn=.true.,fail
2645 real(kind=8),dimension(3) :: e1,e2,e3
2646 real(kind=8) :: dcj,efree_temp
2647 character(len=3) :: seq,res
2648 character(len=5) :: atom
2649 character(len=80) :: card
2650 real(kind=8),dimension(3,20) :: sccor
2651 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
2652 integer :: isugar,molecprev,firstion
2653 character*1 :: sugar
2655 real(kind=8),dimension(3) :: ccc
2657 integer,dimension(2,maxres/3) :: hfrag_alloc
2658 integer,dimension(4,maxres/3) :: bfrag_alloc
2659 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2660 real(kind=8),dimension(:,:), allocatable :: c_temporary
2661 integer,dimension(:,:) , allocatable :: itype_temporary
2662 integer,dimension(:),allocatable :: istype_temp
2669 ! write (2,*) "UNRES_PDB",unres_pdb
2677 !-----------------------------
2678 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2679 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2682 read (ipdbin,'(a80)',end=10) card
2683 ! write (iout,'(a)') card
2684 if (card(:5).eq.'HELIX') then
2687 read(card(22:25),*) hfrag(1,nhfrag)
2688 read(card(34:37),*) hfrag(2,nhfrag)
2690 if (card(:5).eq.'SHEET') then
2693 read(card(24:26),*) bfrag(1,nbfrag)
2694 read(card(35:37),*) bfrag(2,nbfrag)
2695 !rc----------------------------------------
2696 !rc to be corrected !!!
2697 bfrag(3,nbfrag)=bfrag(1,nbfrag)
2698 bfrag(4,nbfrag)=bfrag(2,nbfrag)
2699 !rc----------------------------------------
2701 if (card(:3).eq.'END') then
2703 else if (card(:3).eq.'TER') then
2708 itype(ires_old,molecule)=ntyp1_molec(molecule)
2709 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
2710 nres_molec(molecule)=nres_molec(molecule)+2
2712 ! write (iout,*) "Chain ended",ires,ishift,ires_old
2715 dc(j,ires)=sccor(j,iii)
2718 call sccenter(ires,iii,sccor)
2724 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2725 ! Fish out the ATOM cards.
2726 ! write(iout,*) 'card',card(1:20)
2727 if (index(card(1:4),'ATOM').gt.0) then
2728 read (card(12:16),*) atom
2729 ! write (iout,*) "! ",atom," !",ires
2730 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
2731 read (card(23:26),*) ires
2732 read (card(18:20),'(a3)') res
2733 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
2734 ! & " ires_old",ires_old
2735 ! write (iout,*) "ishift",ishift," ishift1",ishift1
2736 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2737 if (ires-ishift+ishift1.ne.ires_old) then
2738 ! Calculate the CM of the preceding residue.
2739 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2741 ! write (iout,*) "Calculating sidechain center iii",iii
2744 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2747 call sccenter(ires_old,iii,sccor)
2751 ! Start new residue.
2752 if (res.eq.'Cl-' .or. res.eq.'Na+') then
2755 else if (ibeg.eq.1) then
2756 write (iout,*) "BEG ires",ires
2758 if (res.ne.'GLY' .and. res.ne. 'ACE') then
2761 nres_molec(molecule)=nres_molec(molecule)+1
2763 ires=ires-ishift+ishift1
2765 ! write (iout,*) "ishift",ishift," ires",ires,&
2766 ! " ires_old",ires_old
2768 else if (ibeg.eq.2) then
2770 ishift=-ires_old+ires-1 !!!!!
2771 ishift1=ishift1-1 !!!!!
2772 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2773 ires=ires-ishift+ishift1
2774 ! print *,ires,ishift,ishift1
2778 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2779 ires=ires-ishift+ishift1
2782 ! print *,'atom',ires,atom
2783 if (res.eq.'ACE' .or. res.eq.'NHE') then
2786 if (atom.eq.'CA '.or.atom.eq.'N ') then
2788 itype(ires,molecule)=rescode(ires,res,0,molecule)
2790 ! nres_molec(molecule)=nres_molec(molecule)+1
2793 itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
2794 ! nres_molec(molecule)=nres_molec(molecule)+1
2795 read (card(19:19),'(a1)') sugar
2796 isugar=sugarcode(sugar,ires)
2797 ! if (ibeg.eq.1) then
2801 ! print *,"ires=",ires,istype(ires)
2807 ires=ires-ishift+ishift1
2809 ! write (iout,*) "ires_old",ires_old," ires",ires
2810 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2813 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
2814 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2815 res.eq.'NHE'.and.atom(:2).eq.'HN') then
2816 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2817 ! print *,ires,ishift,ishift1
2818 ! write (iout,*) "backbone ",atom
2820 write (iout,'(2i3,2x,a,3f8.3)') &
2821 ires,itype(ires,1),res,(c(j,ires),j=1,3)
2824 nres_molec(molecule)=nres_molec(molecule)+1
2826 sccor(j,iii)=c(j,ires)
2828 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
2829 atom.eq."C2'" .or. atom.eq."C3'" &
2830 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
2831 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
2832 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
2833 ! print *,ires,ishift,ishift1
2837 ! sccor(j,iii)=c(j,ires)
2840 c(j,ires)=c(j,ires)+ccc(j)/5.0
2842 print *,counter,molecule
2843 if (counter.eq.5) then
2845 nres_molec(molecule)=nres_molec(molecule)+1
2848 ! sccor(j,iii)=c(j,ires)
2852 ! print *, "ATOM",atom(1:3)
2853 else if (atom.eq."C5'") then
2854 read (card(19:19),'(a1)') sugar
2855 isugar=sugarcode(sugar,ires)
2860 ! print *,ires,istype(ires)
2863 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2866 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2868 ! write (*,*) card(23:27),ires,itype(ires,1)
2869 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2870 atom.ne.'N' .and. atom.ne.'C' .and. &
2871 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2872 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
2873 .and. atom.ne.'P '.and. &
2874 atom(1:1).ne.'H' .and. &
2875 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
2877 ! write (iout,*) "sidechain ",atom
2878 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2879 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
2880 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2882 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2885 else if ((ions).and.(card(1:6).eq.'HETATM')) then
2886 if (firstion.eq.0) then
2890 dc(j,ires)=sccor(j,iii)
2893 call sccenter(ires,iii,sccor)
2896 read (card(12:16),*) atom
2897 print *,"HETATOM", atom
2898 read (card(18:20),'(a3)') res
2899 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
2900 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
2901 .or.(atom(1:2).eq.'K ')) &
2904 if (molecule.ne.5) molecprev=molecule
2906 nres_molec(molecule)=nres_molec(molecule)+1
2907 print *,"HERE",nres_molec(molecule)
2908 itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
2909 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2913 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
2914 if (ires.eq.0) return
2915 ! Calculate dummy residue coordinates inside the "chain" of a multichain
2918 if ((ires_old.ne.ires).and.(molecule.ne.5)) &
2919 nres_molec(molecule)=nres_molec(molecule)-2
2920 ! print *,'I have', nres_molec(:)
2922 do k=1,4 ! ions are without dummy
2923 if (nres_molec(k).eq.0) cycle
2925 ! write (iout,*) i,itype(i,1)
2926 ! if (itype(i,1).eq.ntyp1) then
2927 ! write (iout,*) "dummy",i,itype(i,1)
2929 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
2930 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
2934 if (itype(i,k).eq.ntyp1_molec(k)) then
2935 if (itype(i+1,k).eq.ntyp1_molec(k)) then
2936 if (itype(i+2,k).eq.0) then
2937 ! print *,"masz sieczke"
2939 if (itype(i+2,j).ne.0) then
2941 itype(i+1,j)=ntyp1_molec(j)
2942 nres_molec(k)=nres_molec(k)-1
2943 nres_molec(j)=nres_molec(j)+1
2949 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
2950 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
2951 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
2953 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2954 ! print *,i,'tu dochodze'
2955 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
2963 c(j,i)=c(j,i-1)-1.9d0*e2(j)
2967 dcj=(c(j,i-2)-c(j,i-3))/2.0
2968 if (dcj.eq.0) dcj=1.23591524223
2973 else !itype(i+1,1).eq.ntyp1
2975 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2976 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
2983 c(j,i)=c(j,i+1)-1.9d0*e2(j)
2987 dcj=(c(j,i+3)-c(j,i+2))/2.0
2988 if (dcj.eq.0) dcj=1.23591524223
2993 endif !itype(i+1,1).eq.ntyp1
2994 endif !itype.eq.ntyp1
2998 ! Calculate the CM of the last side chain.
3002 dc(j,ires)=sccor(j,iii)
3005 call sccenter(ires,iii,sccor)
3011 ! print *,"molecule",molecule
3012 if ((itype(nres,1).ne.10)) then
3014 if (molecule.eq.5) molecule=molecprev
3015 itype(nres,molecule)=ntyp1_molec(molecule)
3016 nres_molec(molecule)=nres_molec(molecule)+1
3018 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3019 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3026 c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3030 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3031 c(j,nres)=c(j,nres-1)+dcj
3032 c(j,2*nres)=c(j,nres)
3036 ! print *,'I have',nres, nres_molec(:)
3038 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3040 if (nres.ne.nres0) then
3041 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3043 stop "Error nres value in WHAM input"
3046 !---------------------------------
3047 !el reallocate tables
3050 ! hfrag_alloc(j,i)=hfrag(j,i)
3053 ! bfrag_alloc(j,i)=bfrag(j,i)
3059 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
3060 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3061 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3062 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
3066 ! hfrag(j,i)=hfrag_alloc(j,i)
3071 ! bfrag(j,i)=bfrag_alloc(j,i)
3074 !el end reallocate tables
3075 !---------------------------------
3083 c(j,2*nres)=c(j,nres)
3086 if (itype(1,1).eq.ntyp1) then
3090 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3091 call refsys(2,3,4,e1,e2,e3,fail)
3098 c(j,1)=c(j,2)-1.9d0*e2(j)
3102 dcj=(c(j,4)-c(j,3))/2.0
3108 ! First lets assign correct dummy to correct type of chain
3110 if (itype(1,1).eq.ntyp1) then
3111 if (itype(2,1).eq.0) then
3113 if (itype(2,j).ne.0) then
3115 itype(1,j)=ntyp1_molec(j)
3116 nres_molec(1)=nres_molec(1)-1
3117 nres_molec(j)=nres_molec(j)+1
3124 print *,'I have',nres, nres_molec(:)
3126 ! Copy the coordinates to reference coordinates
3132 ! Calculate internal coordinates.
3134 write (iout,'(/a)') &
3135 "Cartesian coordinates of the reference structure"
3136 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3137 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3139 write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3140 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3141 (c(j,ires+nres),j=1,3)
3144 ! znamy już nres wiec mozna alokowac tablice
3145 ! Calculate internal coordinates.
3146 if(me.eq.king.or..not.out1file)then
3147 write (iout,'(a)') &
3148 "Backbone and SC coordinates as read from the PDB"
3150 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
3151 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3152 (c(j,nres+ires),j=1,3)
3155 ! NOW LETS ROCK! SORTING
3156 allocate(c_temporary(3,2*nres))
3157 allocate(itype_temporary(nres,5))
3158 allocate(molnum(nres))
3159 allocate(istype_temp(nres))
3160 itype_temporary(:,:)=0
3164 if (itype(i,k).ne.0) then
3166 c_temporary(j,seqalingbegin)=c(j,i)
3167 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3170 itype_temporary(seqalingbegin,k)=itype(i,k)
3171 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
3172 istype_temp(seqalingbegin)=istype(i)
3173 molnum(seqalingbegin)=k
3174 seqalingbegin=seqalingbegin+1
3180 c(j,i)=c_temporary(j,i)
3185 itype(i,k)=itype_temporary(i,k)
3186 istype(i)=istype_temp(i)
3189 if (itype(1,1).eq.ntyp1) then
3193 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3194 call refsys(2,3,4,e1,e2,e3,fail)
3201 c(j,1)=c(j,2)-1.9d0*e2(j)
3205 dcj=(c(j,4)-c(j,3))/2.0
3213 write (iout,'(/a)') &
3214 "Cartesian coordinates of the reference structure after sorting"
3215 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3216 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3218 write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3219 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3220 (c(j,ires+nres),j=1,3)
3224 ! print *,seqalingbegin,nres
3225 if(.not.allocated(vbld)) then
3226 allocate(vbld(2*nres))
3231 if(.not.allocated(vbld_inv)) then
3232 allocate(vbld_inv(2*nres))
3238 if(.not.allocated(theta)) then
3239 allocate(theta(nres+2))
3243 if(.not.allocated(phi)) allocate(phi(nres+2))
3244 if(.not.allocated(alph)) allocate(alph(nres+2))
3245 if(.not.allocated(omeg)) allocate(omeg(nres+2))
3246 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3247 if(.not.allocated(phiref)) allocate(phiref(nres+2))
3248 if(.not.allocated(costtab)) allocate(costtab(nres))
3249 if(.not.allocated(sinttab)) allocate(sinttab(nres))
3250 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3251 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3252 if(.not.allocated(xxref)) allocate(xxref(nres))
3253 if(.not.allocated(yyref)) allocate(yyref(nres))
3254 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3255 if(.not.allocated(dc_norm)) then
3256 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
3257 allocate(dc_norm(3,0:2*nres+2))
3261 call int_from_cart(.true.,.false.)
3262 call sc_loc_geom(.false.)
3264 thetaref(i)=theta(i)
3274 dc(j,i)=c(j,i+1)-c(j,i)
3275 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
3280 dc(j,i+nres)=c(j,i+nres)-c(j,i)
3281 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
3283 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
3287 ! Copy the coordinates to reference coordinates
3288 ! Splits to single chain if occurs
3292 ! cref(j,i,cou)=c(j,i)
3296 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
3297 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
3298 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
3299 !-----------------------------
3303 write (iout,*) "symetr", symetr
3306 !c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3308 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
3311 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3316 cref(j,i,cou)=c(j,i)
3317 cref(j,i+nres,cou)=c(j,i+nres)
3319 chain_rep(j,lll,kkk)=c(j,i)
3320 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
3324 write (iout,*) chain_length
3325 if (chain_length.eq.0) chain_length=nres
3327 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
3328 chain_rep(j,chain_length+nres,symetr) &
3329 =chain_rep(j,chain_length+nres,1)
3332 ! write (iout,*) "spraw lancuchy",chain_length,symetr
3334 ! do kkk=1,chain_length
3335 ! write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
3339 ! makes copy of chains
3340 write (iout,*) "symetr", symetr
3345 if (symetr.gt.1) then
3352 write(iout,*) (tabperm(i,kkk),kkk=1,4)
3358 ! write (iout,*) i,icha
3359 do lll=1,chain_length
3361 if (cou.le.nres) then
3363 kupa=mod(lll,chain_length)
3364 iprzes=(kkk-1)*chain_length+lll
3365 if (kupa.eq.0) kupa=chain_length
3366 ! write (iout,*) "kupa", kupa
3367 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
3368 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
3375 !-koniec robienia kopii
3378 write (iout,*) "nowa struktura", nperm
3380 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
3382 cref(3,i,kkk),cref(1,nres+i,kkk),&
3383 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
3385 100 format (//' alpha-carbon coordinates ',&
3386 ' centroid coordinates'/ &
3387 ' ', 6X,'X',11X,'Y',11X,'Z', &
3388 10X,'X',11X,'Y',11X,'Z')
3389 110 format (a,'(',i3,')',6f12.5)
3395 bfrag(i,j)=bfrag(i,j)-ishift
3401 hfrag(i,j)=hfrag(i,j)-ishift
3407 end subroutine readpdb
3408 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3409 !-----------------------------------------------------------------------------
3411 !-----------------------------------------------------------------------------
3412 subroutine read_control
3426 use random, only: random_init
3427 ! implicit real*8 (a-h,o-z)
3428 ! include 'DIMENSIONS'
3430 use prng, only:prng_restart
3432 logical :: OKRandom!, prng_restart
3435 ! include 'COMMON.IOUNITS'
3436 ! include 'COMMON.TIME1'
3437 ! include 'COMMON.THREAD'
3438 ! include 'COMMON.SBRIDGE'
3439 ! include 'COMMON.CONTROL'
3440 ! include 'COMMON.MCM'
3441 ! include 'COMMON.MAP'
3442 ! include 'COMMON.HEADER'
3443 ! include 'COMMON.CSA'
3444 ! include 'COMMON.CHAIN'
3445 ! include 'COMMON.MUCA'
3446 ! include 'COMMON.MD'
3447 ! include 'COMMON.FFIELD'
3448 ! include 'COMMON.INTERACT'
3449 ! include 'COMMON.SETUP'
3450 !el integer :: KDIAG,ICORFL,IXDR
3451 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
3452 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
3453 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
3454 ! character(len=80) :: ucase
3455 character(len=640) :: controlcard
3457 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
3463 read (INP,'(a)') titel
3464 call card_concat(controlcard,.true.)
3465 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
3466 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
3467 call reada(controlcard,'SEED',seed,0.0D0)
3468 call random_init(seed)
3469 ! Set up the time limit (caution! The time must be input in minutes!)
3470 read_cart=index(controlcard,'READ_CART').gt.0
3471 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3472 call readi(controlcard,'SYM',symetr,1)
3473 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
3474 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
3475 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
3476 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
3477 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
3478 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
3479 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
3480 call reada(controlcard,'DRMS',drms,0.1D0)
3481 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3482 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
3483 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
3484 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
3485 write (iout,'(a,f10.1)')'DRMS = ',drms
3486 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
3487 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
3489 call readi(controlcard,'NZ_START',nz_start,0)
3490 call readi(controlcard,'NZ_END',nz_end,0)
3491 call readi(controlcard,'IZ_SC',iz_sc,0)
3492 timlim=60.0D0*timlim
3493 safety = 60.0d0*safety
3496 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3497 !C SHIELD keyword sets if the shielding effect of side-chains is used
3498 !C 0 denots no shielding is used all peptide are equally despite the
3499 !C solvent accesible area
3500 !C 1 the newly introduced function
3501 !C 2 reseved for further possible developement
3502 call readi(controlcard,'SHIELD',shield_mode,0)
3503 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3504 write(iout,*) "shield_mode",shield_mode
3505 !C Varibles set size of box
3506 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
3507 protein=index(controlcard,"PROTEIN").gt.0
3508 ions=index(controlcard,"IONS").gt.0
3509 nucleic=index(controlcard,"NUCLEIC").gt.0
3510 write (iout,*) "with_theta_constr ",with_theta_constr
3511 AFMlog=(index(controlcard,'AFM'))
3512 selfguide=(index(controlcard,'SELFGUIDE'))
3513 print *,'AFMlog',AFMlog,selfguide,"KUPA"
3514 call readi(controlcard,'GENCONSTR',genconstr,0)
3515 call reada(controlcard,'BOXX',boxxsize,100.0d0)
3516 call reada(controlcard,'BOXY',boxysize,100.0d0)
3517 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
3518 call readi(controlcard,'TUBEMOD',tubemode,0)
3519 write (iout,*) TUBEmode,"TUBEMODE"
3520 if (TUBEmode.gt.0) then
3521 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
3522 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
3523 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
3524 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
3525 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
3526 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
3527 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
3528 buftubebot=bordtubebot+tubebufthick
3529 buftubetop=bordtubetop-tubebufthick
3532 ! CUTOFFF ON ELECTROSTATICS
3533 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
3534 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
3535 write(iout,*) "R_CUT_ELE=",r_cut_ele
3536 ! Lipidic parameters
3537 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3538 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3539 if (lipthick.gt.0.0d0) then
3540 bordliptop=(boxzsize+lipthick)/2.0
3541 bordlipbot=bordliptop-lipthick
3542 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3543 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3544 buflipbot=bordlipbot+lipbufthick
3545 bufliptop=bordliptop-lipbufthick
3546 if ((lipbufthick*2.0d0).gt.lipthick) &
3547 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3548 endif !lipthick.gt.0
3549 write(iout,*) "bordliptop=",bordliptop
3550 write(iout,*) "bordlipbot=",bordlipbot
3551 write(iout,*) "bufliptop=",bufliptop
3552 write(iout,*) "buflipbot=",buflipbot
3553 write (iout,*) "SHIELD MODE",shield_mode
3555 !C-------------------------
3556 minim=(index(controlcard,'MINIMIZE').gt.0)
3557 dccart=(index(controlcard,'CART').gt.0)
3558 overlapsc=(index(controlcard,'OVERLAP').gt.0)
3559 overlapsc=.not.overlapsc
3560 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3561 searchsc=.not.searchsc
3562 sideadd=(index(controlcard,'SIDEADD').gt.0)
3563 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3564 outpdb=(index(controlcard,'PDBOUT').gt.0)
3565 outmol2=(index(controlcard,'MOL2OUT').gt.0)
3566 pdbref=(index(controlcard,'PDBREF').gt.0)
3567 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3568 indpdb=index(controlcard,'PDBSTART')
3569 extconf=(index(controlcard,'EXTCONF').gt.0)
3570 call readi(controlcard,'IPRINT',iprint,0)
3571 call readi(controlcard,'MAXGEN',maxgen,10000)
3572 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3573 call readi(controlcard,"KDIAG",kdiag,0)
3574 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3575 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3576 write (iout,*) "RESCALE_MODE",rescale_mode
3577 split_ene=index(controlcard,'SPLIT_ENE').gt.0
3578 if (index(controlcard,'REGULAR').gt.0.0D0) then
3579 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3583 if (index(controlcard,'CHECKGRAD').gt.0) then
3585 if (index(controlcard,'CART').gt.0) then
3587 elseif (index(controlcard,'CARINT').gt.0) then
3592 elseif (index(controlcard,'THREAD').gt.0) then
3594 call readi(controlcard,'THREAD',nthread,0)
3595 if (nthread.gt.0) then
3596 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3599 write (iout,'(a)')'A number has to follow the THREAD keyword.'
3600 stop 'Error termination in Read_Control.'
3602 else if (index(controlcard,'MCMA').gt.0) then
3604 else if (index(controlcard,'MCEE').gt.0) then
3606 else if (index(controlcard,'MULTCONF').gt.0) then
3608 else if (index(controlcard,'MAP').gt.0) then
3610 call readi(controlcard,'MAP',nmap,0)
3611 else if (index(controlcard,'CSA').gt.0) then
3613 !rc else if (index(controlcard,'ZSCORE').gt.0) then
3615 !rc ZSCORE is rm from UNRES, modecalc=9 is available
3618 !fcm else if (index(controlcard,'MCMF').gt.0) then
3620 else if (index(controlcard,'SOFTREG').gt.0) then
3622 else if (index(controlcard,'CHECK_BOND').gt.0) then
3624 else if (index(controlcard,'TEST').gt.0) then
3626 else if (index(controlcard,'MD').gt.0) then
3628 else if (index(controlcard,'RE ').gt.0) then
3632 lmuca=index(controlcard,'MUCA').gt.0
3633 call readi(controlcard,'MUCADYN',mucadyn,0)
3634 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3635 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3637 write (iout,*) 'MUCADYN=',mucadyn
3638 write (iout,*) 'MUCASMOOTH=',muca_smooth
3641 iscode=index(controlcard,'ONE_LETTER')
3642 indphi=index(controlcard,'PHI')
3643 indback=index(controlcard,'BACK')
3644 iranconf=index(controlcard,'RAND_CONF')
3645 i2ndstr=index(controlcard,'USE_SEC_PRED')
3646 gradout=index(controlcard,'GRADOUT').gt.0
3647 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3648 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3649 if (me.eq.king .or. .not.out1file ) &
3650 write (iout,*) "DISTCHAINMAX",distchainmax
3652 if(me.eq.king.or..not.out1file) &
3653 write (iout,'(2a)') diagmeth(kdiag),&
3654 ' routine used to diagonalize matrices.'
3655 if (shield_mode.gt.0) then
3657 !C VSolvSphere the volume of solving sphere
3659 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3660 !C there will be no distinction between proline peptide group and normal peptide
3661 !C group in case of shielding parameters
3662 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3663 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3664 write (iout,*) VSolvSphere,VSolvSphere_div
3665 !C long axis of side chain
3667 long_r_sidechain(i)=vbldsc0(1,i)
3668 short_r_sidechain(i)=sigma0(i)
3669 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3675 end subroutine read_control
3676 !-----------------------------------------------------------------------------
3677 subroutine read_REMDpar
3679 ! Read REMD settings
3686 use control_data, only:out1file
3687 ! implicit real*8 (a-h,o-z)
3688 ! include 'DIMENSIONS'
3689 ! include 'COMMON.IOUNITS'
3690 ! include 'COMMON.TIME1'
3691 ! include 'COMMON.MD'
3694 !el include 'COMMON.LANGEVIN'
3696 !el include 'COMMON.LANGEVIN.lang0'
3698 ! include 'COMMON.INTERACT'
3699 ! include 'COMMON.NAMES'
3700 ! include 'COMMON.GEO'
3701 ! include 'COMMON.REMD'
3702 ! include 'COMMON.CONTROL'
3703 ! include 'COMMON.SETUP'
3704 ! character(len=80) :: ucase
3705 character(len=320) :: controlcard
3706 character(len=3200) :: controlcard1
3707 integer :: iremd_m_total
3710 ! real(kind=8) :: var,ene
3712 if(me.eq.king.or..not.out1file) &
3713 write (iout,*) "REMD setup"
3715 call card_concat(controlcard,.true.)
3716 call readi(controlcard,"NREP",nrep,3)
3717 call readi(controlcard,"NSTEX",nstex,1000)
3718 call reada(controlcard,"RETMIN",retmin,10.0d0)
3719 call reada(controlcard,"RETMAX",retmax,1000.0d0)
3720 mremdsync=(index(controlcard,'SYNC').gt.0)
3721 call readi(controlcard,"NSYN",i_sync_step,100)
3722 restart1file=(index(controlcard,'REST1FILE').gt.0)
3723 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3724 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3725 if(max_cache_traj_use.gt.max_cache_traj) &
3726 max_cache_traj_use=max_cache_traj
3727 if(me.eq.king.or..not.out1file) then
3728 !d if (traj1file) then
3729 !rc caching is in testing - NTWX is not ignored
3730 !d write (iout,*) "NTWX value is ignored"
3731 !d write (iout,*) " trajectory is stored to one file by master"
3732 !d write (iout,*) " before exchange at NSTEX intervals"
3734 write (iout,*) "NREP= ",nrep
3735 write (iout,*) "NSTEX= ",nstex
3736 write (iout,*) "SYNC= ",mremdsync
3737 write (iout,*) "NSYN= ",i_sync_step
3738 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3741 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3742 if (index(controlcard,'TLIST').gt.0) then
3744 call card_concat(controlcard1,.true.)
3745 read(controlcard1,*) (remd_t(i),i=1,nrep)
3746 if(me.eq.king.or..not.out1file) &
3747 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3750 if (index(controlcard,'MLIST').gt.0) then
3752 call card_concat(controlcard1,.true.)
3753 read(controlcard1,*) (remd_m(i),i=1,nrep)
3754 if(me.eq.king.or..not.out1file) then
3755 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3758 iremd_m_total=iremd_m_total+remd_m(i)
3760 write (iout,*) 'Total number of replicas ',iremd_m_total
3763 if(me.eq.king.or..not.out1file) &
3764 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3766 end subroutine read_REMDpar
3767 !-----------------------------------------------------------------------------
3768 subroutine read_MDpar
3772 use control_data, only: r_cut,rlamb,out1file
3774 use geometry_data, only: pi
3776 ! implicit real*8 (a-h,o-z)
3777 ! include 'DIMENSIONS'
3778 ! include 'COMMON.IOUNITS'
3779 ! include 'COMMON.TIME1'
3780 ! include 'COMMON.MD'
3783 !el include 'COMMON.LANGEVIN'
3785 !el include 'COMMON.LANGEVIN.lang0'
3787 ! include 'COMMON.INTERACT'
3788 ! include 'COMMON.NAMES'
3789 ! include 'COMMON.GEO'
3790 ! include 'COMMON.SETUP'
3791 ! include 'COMMON.CONTROL'
3792 ! include 'COMMON.SPLITELE'
3793 ! character(len=80) :: ucase
3794 character(len=320) :: controlcard
3799 call card_concat(controlcard,.true.)
3800 call readi(controlcard,"NSTEP",n_timestep,1000000)
3801 call readi(controlcard,"NTWE",ntwe,100)
3802 call readi(controlcard,"NTWX",ntwx,1000)
3803 call reada(controlcard,"DT",d_time,1.0d-1)
3804 call reada(controlcard,"DVMAX",dvmax,2.0d1)
3805 call reada(controlcard,"DAMAX",damax,1.0d1)
3806 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3807 call readi(controlcard,"LANG",lang,0)
3808 RESPA = index(controlcard,"RESPA") .gt. 0
3809 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3810 ntime_split0=ntime_split
3811 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3812 ntime_split0=ntime_split
3813 call reada(controlcard,"R_CUT",r_cut,2.0d0)
3814 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3815 rest = index(controlcard,"REST").gt.0
3816 tbf = index(controlcard,"TBF").gt.0
3817 usampl = index(controlcard,"USAMPL").gt.0
3818 mdpdb = index(controlcard,"MDPDB").gt.0
3819 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3820 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3821 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3822 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3823 if (count_reset_moment.eq.0) count_reset_moment=1000000000
3824 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3825 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3826 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3827 if (count_reset_vel.eq.0) count_reset_vel=1000000000
3828 large = index(controlcard,"LARGE").gt.0
3829 print_compon = index(controlcard,"PRINT_COMPON").gt.0
3830 rattle = index(controlcard,"RATTLE").gt.0
3831 ! if performing umbrella sampling, fragments constrained are read from the fragment file
3837 if(me.eq.king.or..not.out1file) then
3839 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3841 write (iout,'(a)') "The units are:"
3842 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3843 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3844 " acceleration: angstrom/(48.9 fs)**2"
3845 write (iout,'(a)') "energy: kcal/mol, temperature: K"
3847 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3848 write (iout,'(a60,f10.5,a)') &
3849 "Initial time step of numerical integration:",d_time,&
3851 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3853 write (iout,'(2a,i4,a)') &
3854 "A-MTS algorithm used; initial time step for fast-varying",&
3855 " short-range forces split into",ntime_split," steps."
3856 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3857 r_cut," lambda",rlamb
3859 write (iout,'(2a,f10.5)') &
3860 "Maximum acceleration threshold to reduce the time step",&
3861 "/increase split number:",damax
3862 write (iout,'(2a,f10.5)') &
3863 "Maximum predicted energy drift to reduce the timestep",&
3864 "/increase split number:",edriftmax
3865 write (iout,'(a60,f10.5)') &
3866 "Maximum velocity threshold to reduce velocities:",dvmax
3867 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3868 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3869 if (rattle) write (iout,'(a60)') &
3870 "Rattle algorithm used to constrain the virtual bonds"
3874 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3875 call reada(controlcard,"RWAT",rwat,1.4d0)
3876 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3877 surfarea=index(controlcard,"SURFAREA").gt.0
3878 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3879 if(me.eq.king.or..not.out1file)then
3880 write (iout,'(/a,$)') "Langevin dynamics calculation"
3882 write (iout,'(a/)') &
3883 " with direct integration of Langevin equations"
3884 else if (lang.eq.2) then
3885 write (iout,'(a/)') " with TINKER stochasic MD integrator"
3886 else if (lang.eq.3) then
3887 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3888 else if (lang.eq.4) then
3889 write (iout,'(a/)') " in overdamped mode"
3891 write (iout,'(//a,i5)') &
3892 "=========== ERROR: Unknown Langevin dynamics mode:",lang
3895 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3896 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
3897 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
3898 write (iout,'(a60,f10.5)') &
3899 "Scaling factor of the friction forces:",scal_fric
3900 if (surfarea) write (iout,'(2a,i10,a)') &
3901 "Friction coefficients will be scaled by solvent-accessible",&
3902 " surface area every",reset_fricmat," steps."
3904 ! Calculate friction coefficients and bounds of stochastic forces
3905 eta=6*pi*cPoise*etawat
3906 if(me.eq.king.or..not.out1file) &
3907 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
3910 do j=1,5 !types of molecules
3911 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
3912 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
3914 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
3915 do j=1,5 !types of molecules
3917 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
3918 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
3922 if(me.eq.king.or..not.out1file)then
3923 write (iout,'(/2a/)') &
3924 "Radii of site types and friction coefficients and std's of",&
3925 " stochastic forces of fully exposed sites"
3926 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
3928 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
3929 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
3933 if(me.eq.king.or..not.out1file)then
3934 write (iout,'(a)') "Berendsen bath calculation"
3935 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3936 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
3938 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
3939 count_reset_moment," steps"
3941 write (iout,'(a,i10,a)') &
3942 "Velocities will be reset at random every",count_reset_vel,&
3946 if(me.eq.king.or..not.out1file) &
3947 write (iout,'(a31)') "Microcanonical mode calculation"
3949 if(me.eq.king.or..not.out1file)then
3950 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
3952 write(iout,*) "MD running with constraints."
3953 write(iout,*) "Equilibration time ", eq_time, " mtus."
3954 write(iout,*) "Constraining ", nfrag," fragments."
3955 write(iout,*) "Length of each fragment, weight and q0:"
3957 write (iout,*) "Set of restraints #",iset
3959 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
3960 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
3962 write(iout,*) "constraints between ", npair, "fragments."
3963 write(iout,*) "constraint pairs, weights and q0:"
3965 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
3966 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
3968 write(iout,*) "angle constraints within ", nfrag_back,&
3969 "backbone fragments."
3970 write(iout,*) "fragment, weights:"
3972 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
3973 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
3974 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
3977 iset=mod(kolor,nset)+1
3980 if(me.eq.king.or..not.out1file) &
3981 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
3983 end subroutine read_MDpar
3984 !-----------------------------------------------------------------------------
3988 ! implicit real*8 (a-h,o-z)
3989 ! include 'DIMENSIONS'
3990 ! include 'COMMON.MAP'
3991 ! include 'COMMON.IOUNITS'
3992 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
3993 character(len=80) :: mapcard !,ucase
3996 ! real(kind=8) :: var,ene
3999 read (inp,'(a)') mapcard
4000 mapcard=ucase(mapcard)
4001 if (index(mapcard,'PHI').gt.0) then
4003 else if (index(mapcard,'THE').gt.0) then
4005 else if (index(mapcard,'ALP').gt.0) then
4007 else if (index(mapcard,'OME').gt.0) then
4010 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4011 stop 'Error - illegal variable spec in MAP card.'
4013 call readi (mapcard,'RES1',res1(imap),0)
4014 call readi (mapcard,'RES2',res2(imap),0)
4015 if (res1(imap).eq.0) then
4016 res1(imap)=res2(imap)
4017 else if (res2(imap).eq.0) then
4018 res2(imap)=res1(imap)
4020 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4021 write (iout,'(a)') &
4022 'Error - illegal definition of variable group in MAP.'
4023 stop 'Error - illegal definition of variable group in MAP.'
4025 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4026 call reada(mapcard,'TO',ang_to(imap),0.0D0)
4027 call readi(mapcard,'NSTEP',nstep(imap),0)
4028 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4029 write (iout,'(a)') &
4030 'Illegal boundary and/or step size specification in MAP.'
4031 stop 'Illegal boundary and/or step size specification in MAP.'
4035 end subroutine map_read
4036 !-----------------------------------------------------------------------------
4039 use control_data, only: vdisulf
4041 ! implicit real*8 (a-h,o-z)
4042 ! include 'DIMENSIONS'
4043 ! include 'COMMON.IOUNITS'
4044 ! include 'COMMON.GEO'
4045 ! include 'COMMON.CSA'
4046 ! include 'COMMON.BANK'
4047 ! include 'COMMON.CONTROL'
4048 ! character(len=80) :: ucase
4049 character(len=620) :: mcmcard
4051 ! integer :: ntf,ik,iw_pdb
4052 ! real(kind=8) :: var,ene
4054 call card_concat(mcmcard,.true.)
4056 call readi(mcmcard,'NCONF',nconf,50)
4057 call readi(mcmcard,'NADD',nadd,0)
4058 call readi(mcmcard,'JSTART',jstart,1)
4059 call readi(mcmcard,'JEND',jend,1)
4060 call readi(mcmcard,'NSTMAX',nstmax,500000)
4061 call readi(mcmcard,'N0',n0,1)
4062 call readi(mcmcard,'N1',n1,6)
4063 call readi(mcmcard,'N2',n2,4)
4064 call readi(mcmcard,'N3',n3,0)
4065 call readi(mcmcard,'N4',n4,0)
4066 call readi(mcmcard,'N5',n5,0)
4067 call readi(mcmcard,'N6',n6,10)
4068 call readi(mcmcard,'N7',n7,0)
4069 call readi(mcmcard,'N8',n8,0)
4070 call readi(mcmcard,'N9',n9,0)
4071 call readi(mcmcard,'N14',n14,0)
4072 call readi(mcmcard,'N15',n15,0)
4073 call readi(mcmcard,'N16',n16,0)
4074 call readi(mcmcard,'N17',n17,0)
4075 call readi(mcmcard,'N18',n18,0)
4077 vdisulf=(index(mcmcard,'DYNSS').gt.0)
4079 call readi(mcmcard,'NDIFF',ndiff,2)
4080 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
4081 call readi(mcmcard,'IS1',is1,1)
4082 call readi(mcmcard,'IS2',is2,8)
4083 call readi(mcmcard,'NRAN0',nran0,4)
4084 call readi(mcmcard,'NRAN1',nran1,2)
4085 call readi(mcmcard,'IRR',irr,1)
4086 call readi(mcmcard,'NSEED',nseed,20)
4087 call readi(mcmcard,'NTOTAL',ntotal,10000)
4088 call reada(mcmcard,'CUT1',cut1,2.0d0)
4089 call reada(mcmcard,'CUT2',cut2,5.0d0)
4090 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
4091 call readi(mcmcard,'ICMAX',icmax,3)
4092 call readi(mcmcard,'IRESTART',irestart,0)
4093 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
4096 call reada(mcmcard,'DELE',dele,20.0d0)
4097 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4098 call readi(mcmcard,'IREF',iref,0)
4099 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4100 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4101 call readi(mcmcard,'NCONF_IN',nconf_in,0)
4102 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4103 write (iout,*) "NCONF_IN",nconf_in
4105 end subroutine csaread
4106 !-----------------------------------------------------------------------------
4110 use control_data, only: MaxMoveType
4113 ! implicit real*8 (a-h,o-z)
4114 ! include 'DIMENSIONS'
4115 ! include 'COMMON.MCM'
4116 ! include 'COMMON.MCE'
4117 ! include 'COMMON.IOUNITS'
4118 ! character(len=80) :: ucase
4119 character(len=320) :: mcmcard
4122 ! real(kind=8) :: var,ene
4124 call card_concat(mcmcard,.true.)
4125 call readi(mcmcard,'MAXACC',maxacc,100)
4126 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4127 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4128 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4129 call readi(mcmcard,'MAXREPM',maxrepm,200)
4130 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4131 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4132 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4133 call reada(mcmcard,'E_UP',e_up,5.0D0)
4134 call reada(mcmcard,'DELTE',delte,0.1D0)
4135 call readi(mcmcard,'NSWEEP',nsweep,5)
4136 call readi(mcmcard,'NSTEPH',nsteph,0)
4137 call readi(mcmcard,'NSTEPC',nstepc,0)
4138 call reada(mcmcard,'TMIN',tmin,298.0D0)
4139 call reada(mcmcard,'TMAX',tmax,298.0D0)
4140 call readi(mcmcard,'NWINDOW',nwindow,0)
4141 call readi(mcmcard,'PRINT_MC',print_mc,0)
4142 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4143 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4144 ent_read=(index(mcmcard,'ENT_READ').gt.0)
4145 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4146 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4147 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4148 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4149 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4150 if (nwindow.gt.0) then
4151 allocate(winstart(nwindow)) !!el (maxres)
4152 allocate(winend(nwindow)) !!el
4153 allocate(winlen(nwindow)) !!el
4154 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4156 winlen(i)=winend(i)-winstart(i)+1
4159 if (tmax.lt.tmin) tmax=tmin
4160 if (tmax.eq.tmin) then
4164 if (nstepc.gt.0 .and. nsteph.gt.0) then
4165 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4166 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4168 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4169 ! Probabilities of different move types
4170 sumpro_type(0)=0.0D0
4171 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4172 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4173 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4174 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
4175 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4176 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4177 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4179 print *,'i',i,' sumprotype',sumpro_type(i)
4180 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4181 print *,'i',i,' sumprotype',sumpro_type(i)
4184 end subroutine mcmread
4185 !-----------------------------------------------------------------------------
4186 subroutine read_minim
4189 ! implicit real*8 (a-h,o-z)
4190 ! include 'DIMENSIONS'
4191 ! include 'COMMON.MINIM'
4192 ! include 'COMMON.IOUNITS'
4193 ! character(len=80) :: ucase
4194 character(len=320) :: minimcard
4196 ! integer :: ntf,ik,iw_pdb
4197 ! real(kind=8) :: var,ene
4199 call card_concat(minimcard,.true.)
4200 call readi(minimcard,'MAXMIN',maxmin,2000)
4201 call readi(minimcard,'MAXFUN',maxfun,5000)
4202 call readi(minimcard,'MINMIN',minmin,maxmin)
4203 call readi(minimcard,'MINFUN',minfun,maxmin)
4204 call reada(minimcard,'TOLF',tolf,1.0D-2)
4205 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4206 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4207 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4208 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4209 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4210 'Options in energy minimization:'
4211 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4212 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4213 'MinMin:',MinMin,' MinFun:',MinFun,&
4214 ' TolF:',TolF,' RTolF:',RTolF
4216 end subroutine read_minim
4217 !-----------------------------------------------------------------------------
4218 subroutine openunits
4220 use MD_data, only: usampl
4223 use control_data, only:out1file
4224 use control, only: getenv_loc
4225 ! implicit real*8 (a-h,o-z)
4226 ! include 'DIMENSIONS'
4229 character(len=16) :: form,nodename
4230 integer :: nodelen,ierror,npos
4232 ! include 'COMMON.SETUP'
4233 ! include 'COMMON.IOUNITS'
4234 ! include 'COMMON.MD'
4235 ! include 'COMMON.CONTROL'
4236 integer :: lenpre,lenpot,lentmp !,ilen
4238 character(len=3) :: out1file_text !,ucase
4239 character(len=3) :: ll
4242 ! integer :: ntf,ik,iw_pdb
4243 ! real(kind=8) :: var,ene
4245 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
4246 call getenv_loc("PREFIX",prefix)
4248 call getenv_loc("POT",pot)
4249 call getenv_loc("DIRTMP",tmpdir)
4250 call getenv_loc("CURDIR",curdir)
4251 call getenv_loc("OUT1FILE",out1file_text)
4252 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
4253 out1file_text=ucase(out1file_text)
4254 if (out1file_text(1:1).eq."Y") then
4257 out1file=fg_rank.gt.0
4262 if (lentmp.gt.0) then
4263 write (*,'(80(1h!))')
4264 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
4265 write (*,'(80(1h!))')
4266 write (*,*)"All output files will be on node /tmp directory."
4268 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
4269 if (me.eq.king) then
4270 write (*,*) "The master node is ",nodename
4271 else if (fg_rank.eq.0) then
4272 write (*,*) "I am the CG slave node ",nodename
4274 write (*,*) "I am the FG slave node ",nodename
4277 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
4278 lenpre = lentmp+lenpre+1
4280 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
4281 ! Get the names and open the input files
4282 #if defined(WINIFL) || defined(WINPGI)
4283 open(1,file=pref_orig(:ilen(pref_orig))// &
4284 '.inp',status='old',readonly,shared)
4285 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4286 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4287 ! Get parameter filenames and open the parameter files.
4288 call getenv_loc('BONDPAR',bondname)
4289 open (ibond,file=bondname,status='old',readonly,shared)
4290 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4291 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
4292 call getenv_loc('THETPAR',thetname)
4293 open (ithep,file=thetname,status='old',readonly,shared)
4294 call getenv_loc('ROTPAR',rotname)
4295 open (irotam,file=rotname,status='old',readonly,shared)
4296 call getenv_loc('TORPAR',torname)
4297 open (itorp,file=torname,status='old',readonly,shared)
4298 call getenv_loc('TORDPAR',tordname)
4299 open (itordp,file=tordname,status='old',readonly,shared)
4300 call getenv_loc('FOURIER',fouriername)
4301 open (ifourier,file=fouriername,status='old',readonly,shared)
4302 call getenv_loc('ELEPAR',elename)
4303 open (ielep,file=elename,status='old',readonly,shared)
4304 call getenv_loc('SIDEPAR',sidename)
4305 open (isidep,file=sidename,status='old',readonly,shared)
4307 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4308 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
4309 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4310 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
4311 call getenv_loc('TORPAR_NUCL',torname_nucl)
4312 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
4313 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4314 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
4315 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4316 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
4319 #elif (defined CRAY) || (defined AIX)
4320 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4322 ! print *,"Processor",myrank," opened file 1"
4323 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4324 ! print *,"Processor",myrank," opened file 9"
4325 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4326 ! Get parameter filenames and open the parameter files.
4327 call getenv_loc('BONDPAR',bondname)
4328 open (ibond,file=bondname,status='old',action='read')
4329 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4330 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4332 ! print *,"Processor",myrank," opened file IBOND"
4333 call getenv_loc('THETPAR',thetname)
4334 open (ithep,file=thetname,status='old',action='read')
4335 ! print *,"Processor",myrank," opened file ITHEP"
4336 call getenv_loc('ROTPAR',rotname)
4337 open (irotam,file=rotname,status='old',action='read')
4338 ! print *,"Processor",myrank," opened file IROTAM"
4339 call getenv_loc('TORPAR',torname)
4340 open (itorp,file=torname,status='old',action='read')
4341 ! print *,"Processor",myrank," opened file ITORP"
4342 call getenv_loc('TORDPAR',tordname)
4343 open (itordp,file=tordname,status='old',action='read')
4344 ! print *,"Processor",myrank," opened file ITORDP"
4345 call getenv_loc('SCCORPAR',sccorname)
4346 open (isccor,file=sccorname,status='old',action='read')
4347 ! print *,"Processor",myrank," opened file ISCCOR"
4348 call getenv_loc('FOURIER',fouriername)
4349 open (ifourier,file=fouriername,status='old',action='read')
4350 ! print *,"Processor",myrank," opened file IFOURIER"
4351 call getenv_loc('ELEPAR',elename)
4352 open (ielep,file=elename,status='old',action='read')
4353 ! print *,"Processor",myrank," opened file IELEP"
4354 call getenv_loc('SIDEPAR',sidename)
4355 open (isidep,file=sidename,status='old',action='read')
4357 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4358 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4359 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4360 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4361 call getenv_loc('TORPAR_NUCL',torname_nucl)
4362 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4363 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4364 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4365 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4366 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4368 call getenv_loc('LIPTRANPAR',liptranname)
4369 open (iliptranpar,file=liptranname,status='old',action='read')
4370 call getenv_loc('TUBEPAR',tubename)
4371 open (itube,file=tubename,status='old',action='read')
4372 call getenv_loc('IONPAR',ionname)
4373 open (iion,file=ionname,status='old',action='read')
4375 ! print *,"Processor",myrank," opened file ISIDEP"
4376 ! print *,"Processor",myrank," opened parameter files"
4378 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
4379 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4380 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4381 ! Get parameter filenames and open the parameter files.
4382 call getenv_loc('BONDPAR',bondname)
4383 open (ibond,file=bondname,status='old')
4384 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4385 open (ibond_nucl,file=bondname_nucl,status='old')
4387 call getenv_loc('THETPAR',thetname)
4388 open (ithep,file=thetname,status='old')
4389 call getenv_loc('ROTPAR',rotname)
4390 open (irotam,file=rotname,status='old')
4391 call getenv_loc('TORPAR',torname)
4392 open (itorp,file=torname,status='old')
4393 call getenv_loc('TORDPAR',tordname)
4394 open (itordp,file=tordname,status='old')
4395 call getenv_loc('SCCORPAR',sccorname)
4396 open (isccor,file=sccorname,status='old')
4397 call getenv_loc('FOURIER',fouriername)
4398 open (ifourier,file=fouriername,status='old')
4399 call getenv_loc('ELEPAR',elename)
4400 open (ielep,file=elename,status='old')
4401 call getenv_loc('SIDEPAR',sidename)
4402 open (isidep,file=sidename,status='old')
4404 open (ithep_nucl,file=thetname_nucl,status='old')
4405 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4406 open (irotam_nucl,file=rotname_nucl,status='old')
4407 call getenv_loc('TORPAR_NUCL',torname_nucl)
4408 open (itorp_nucl,file=torname_nucl,status='old')
4409 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4410 open (itordp_nucl,file=tordname_nucl,status='old')
4411 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4412 open (isidep_nucl,file=sidename_nucl,status='old')
4414 call getenv_loc('LIPTRANPAR',liptranname)
4415 open (iliptranpar,file=liptranname,status='old')
4416 call getenv_loc('TUBEPAR',tubename)
4417 open (itube,file=tubename,status='old')
4418 call getenv_loc('IONPAR',ionname)
4419 open (iion,file=ionname,status='old')
4421 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4423 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4424 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4425 ! Get parameter filenames and open the parameter files.
4426 call getenv_loc('BONDPAR',bondname)
4427 open (ibond,file=bondname,status='old',action='read')
4428 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4429 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4430 call getenv_loc('THETPAR',thetname)
4431 open (ithep,file=thetname,status='old',action='read')
4432 call getenv_loc('ROTPAR',rotname)
4433 open (irotam,file=rotname,status='old',action='read')
4434 call getenv_loc('TORPAR',torname)
4435 open (itorp,file=torname,status='old',action='read')
4436 call getenv_loc('TORDPAR',tordname)
4437 open (itordp,file=tordname,status='old',action='read')
4438 call getenv_loc('SCCORPAR',sccorname)
4439 open (isccor,file=sccorname,status='old',action='read')
4441 call getenv_loc('THETPARPDB',thetname_pdb)
4442 print *,"thetname_pdb ",thetname_pdb
4443 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
4444 print *,ithep_pdb," opened"
4446 call getenv_loc('FOURIER',fouriername)
4447 open (ifourier,file=fouriername,status='old',readonly)
4448 call getenv_loc('ELEPAR',elename)
4449 open (ielep,file=elename,status='old',readonly)
4450 call getenv_loc('SIDEPAR',sidename)
4451 open (isidep,file=sidename,status='old',readonly)
4453 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4454 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4455 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4456 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4457 call getenv_loc('TORPAR_NUCL',torname_nucl)
4458 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4459 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4460 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4461 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4462 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4463 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
4464 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
4466 call getenv_loc('LIPTRANPAR',liptranname)
4467 open (iliptranpar,file=liptranname,status='old',action='read')
4468 call getenv_loc('TUBEPAR',tubename)
4469 open (itube,file=tubename,status='old',action='read')
4470 call getenv_loc('IONPAR',ionname)
4471 open (iion,file=ionname,status='old',action='read')
4474 call getenv_loc('ROTPARPDB',rotname_pdb)
4475 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
4478 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
4479 #if defined(WINIFL) || defined(WINPGI)
4480 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
4481 #elif (defined CRAY) || (defined AIX)
4482 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4484 open (iscpp_nucl,file=scpname_nucl,status='old')
4486 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4491 ! 8/9/01 In the newest version SCp interaction constants are read from a file
4492 ! Use -DOLDSCP to use hard-coded constants instead.
4494 call getenv_loc('SCPPAR',scpname)
4495 #if defined(WINIFL) || defined(WINPGI)
4496 open (iscpp,file=scpname,status='old',readonly,shared)
4497 #elif (defined CRAY) || (defined AIX)
4498 open (iscpp,file=scpname,status='old',action='read')
4500 open (iscpp,file=scpname,status='old')
4502 open (iscpp,file=scpname,status='old',action='read')
4505 call getenv_loc('PATTERN',patname)
4506 #if defined(WINIFL) || defined(WINPGI)
4507 open (icbase,file=patname,status='old',readonly,shared)
4508 #elif (defined CRAY) || (defined AIX)
4509 open (icbase,file=patname,status='old',action='read')
4511 open (icbase,file=patname,status='old')
4513 open (icbase,file=patname,status='old',action='read')
4516 ! Open output file only for CG processes
4517 ! print *,"Processor",myrank," fg_rank",fg_rank
4518 if (fg_rank.eq.0) then
4520 if (nodes.eq.1) then
4523 npos = dlog10(dfloat(nodes-1))+1
4525 if (npos.lt.3) npos=3
4526 write (liczba,'(i1)') npos
4527 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
4529 write (liczba,form) me
4530 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
4531 liczba(:ilen(liczba))
4532 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4534 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4536 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
4537 liczba(:ilen(liczba))//'.mol2'
4538 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4539 liczba(:ilen(liczba))//'.stat'
4541 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
4542 //liczba(:ilen(liczba))//'.stat')
4543 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
4546 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4547 liczba(:ilen(liczba))//'.const'
4552 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
4553 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
4554 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
4555 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
4556 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
4558 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
4560 rest2name=prefix(:ilen(prefix))//'.rst'
4562 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
4565 #if defined(AIX) || defined(PGI)
4566 if (me.eq.king .or. .not. out1file) &
4567 open(iout,file=outname,status='unknown')
4569 if (fg_rank.gt.0) then
4570 write (liczba,'(i3.3)') myrank/nfgtasks
4571 write (ll,'(bz,i3.3)') fg_rank
4572 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4577 open(igeom,file=intname,status='unknown',position='append')
4578 open(ipdb,file=pdbname,status='unknown')
4579 open(imol2,file=mol2name,status='unknown')
4580 open(istat,file=statname,status='unknown',position='append')
4582 !1out open(iout,file=outname,status='unknown')
4585 if (me.eq.king .or. .not.out1file) &
4586 open(iout,file=outname,status='unknown')
4588 if (fg_rank.gt.0) then
4589 write (liczba,'(i3.3)') myrank/nfgtasks
4590 write (ll,'(bz,i3.3)') fg_rank
4591 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4596 open(igeom,file=intname,status='unknown',access='append')
4597 open(ipdb,file=pdbname,status='unknown')
4598 open(imol2,file=mol2name,status='unknown')
4599 open(istat,file=statname,status='unknown',access='append')
4601 !1out open(iout,file=outname,status='unknown')
4604 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
4605 csa_seed=prefix(:lenpre)//'.CSA.seed'
4606 csa_history=prefix(:lenpre)//'.CSA.history'
4607 csa_bank=prefix(:lenpre)//'.CSA.bank'
4608 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
4609 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
4610 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
4611 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
4612 csa_int=prefix(:lenpre)//'.int'
4613 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
4614 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
4615 csa_in=prefix(:lenpre)//'.CSA.in'
4616 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
4619 write (iout,'(80(1h-))')
4620 write (iout,'(30x,a)') "FILE ASSIGNMENT"
4621 write (iout,'(80(1h-))')
4622 write (iout,*) "Input file : ",&
4623 pref_orig(:ilen(pref_orig))//'.inp'
4624 write (iout,*) "Output file : ",&
4625 outname(:ilen(outname))
4627 write (iout,*) "Sidechain potential file : ",&
4628 sidename(:ilen(sidename))
4630 write (iout,*) "SCp potential file : ",&
4631 scpname(:ilen(scpname))
4633 write (iout,*) "Electrostatic potential file : ",&
4634 elename(:ilen(elename))
4635 write (iout,*) "Cumulant coefficient file : ",&
4636 fouriername(:ilen(fouriername))
4637 write (iout,*) "Torsional parameter file : ",&
4638 torname(:ilen(torname))
4639 write (iout,*) "Double torsional parameter file : ",&
4640 tordname(:ilen(tordname))
4641 write (iout,*) "SCCOR parameter file : ",&
4642 sccorname(:ilen(sccorname))
4643 write (iout,*) "Bond & inertia constant file : ",&
4644 bondname(:ilen(bondname))
4645 write (iout,*) "Bending parameter file : ",&
4646 thetname(:ilen(thetname))
4647 write (iout,*) "Rotamer parameter file : ",&
4648 rotname(:ilen(rotname))
4651 write (iout,*) "Thetpdb parameter file : ",&
4652 thetname_pdb(:ilen(thetname_pdb))
4655 write (iout,*) "Threading database : ",&
4656 patname(:ilen(patname))
4658 write (iout,*)" DIRTMP : ",&
4660 write (iout,'(80(1h-))')
4663 end subroutine openunits
4664 !-----------------------------------------------------------------------------
4667 use geometry_data, only: nres,dc
4669 ! implicit real*8 (a-h,o-z)
4670 ! include 'DIMENSIONS'
4671 ! include 'COMMON.CHAIN'
4672 ! include 'COMMON.IOUNITS'
4673 ! include 'COMMON.MD'
4676 ! real(kind=8) :: var,ene
4678 open(irest2,file=rest2name,status='unknown')
4679 read(irest2,*) totT,EK,potE,totE,t_bath
4682 ! AL 4/17/17: Now reading d_t(0,:) too
4684 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4687 ! AL 4/17/17: Now reading d_c(0,:) too
4689 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4692 read (irest2,*) iset
4696 end subroutine readrst
4697 !-----------------------------------------------------------------------------
4698 subroutine read_fragments
4702 use control_data, only:out1file
4705 ! implicit real*8 (a-h,o-z)
4706 ! include 'DIMENSIONS'
4710 ! include 'COMMON.SETUP'
4711 ! include 'COMMON.CHAIN'
4712 ! include 'COMMON.IOUNITS'
4713 ! include 'COMMON.MD'
4714 ! include 'COMMON.CONTROL'
4717 ! real(kind=8) :: var,ene
4719 read(inp,*) nset,nfrag,npair,nfrag_back
4721 !el from module energy
4722 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
4723 if(.not.allocated(wfrag_back)) then
4724 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4725 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4727 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4728 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
4730 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4731 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
4734 if(me.eq.king.or..not.out1file) &
4735 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4736 " nfrag_back",nfrag_back
4738 read(inp,*) mset(iset)
4740 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4742 if(me.eq.king.or..not.out1file) &
4743 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4744 ifrag(2,i,iset), qinfrag(i,iset)
4747 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4749 if(me.eq.king.or..not.out1file) &
4750 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4751 ipair(2,i,iset), qinpair(i,iset)
4754 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4755 wfrag_back(3,i,iset),&
4756 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4757 if(me.eq.king.or..not.out1file) &
4758 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4759 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4763 end subroutine read_fragments
4764 !-----------------------------------------------------------------------------
4766 !-----------------------------------------------------------------------------
4770 ! implicit real*8 (a-h,o-z)
4771 ! include 'DIMENSIONS'
4772 ! include 'COMMON.CSA'
4773 ! include 'COMMON.BANK'
4774 ! include 'COMMON.IOUNITS'
4776 ! integer :: ntf,ik,iw_pdb
4777 ! real(kind=8) :: var,ene
4779 open(icsa_in,file=csa_in,status="old",err=100)
4780 read(icsa_in,*) nconf
4781 read(icsa_in,*) jstart,jend
4782 read(icsa_in,*) nstmax
4783 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4784 read(icsa_in,*) nran0,nran1,irr
4785 read(icsa_in,*) nseed
4786 read(icsa_in,*) ntotal,cut1,cut2
4787 read(icsa_in,*) estop
4788 read(icsa_in,*) icmax,irestart
4789 read(icsa_in,*) ntbankm,dele,difcut
4790 read(icsa_in,*) iref,rmscut,pnccut
4791 read(icsa_in,*) ndiff
4798 end subroutine csa_read
4799 !-----------------------------------------------------------------------------
4800 subroutine initial_write
4803 ! implicit real*8 (a-h,o-z)
4804 ! include 'DIMENSIONS'
4805 ! include 'COMMON.CSA'
4806 ! include 'COMMON.BANK'
4807 ! include 'COMMON.IOUNITS'
4809 ! integer :: ntf,ik,iw_pdb
4810 ! real(kind=8) :: var,ene
4812 open(icsa_seed,file=csa_seed,status="unknown")
4813 write(icsa_seed,*) "seed"
4815 #if defined(AIX) || defined(PGI)
4816 open(icsa_history,file=csa_history,status="unknown",&
4819 open(icsa_history,file=csa_history,status="unknown",&
4822 write(icsa_history,*) nconf
4823 write(icsa_history,*) jstart,jend
4824 write(icsa_history,*) nstmax
4825 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4826 write(icsa_history,*) nran0,nran1,irr
4827 write(icsa_history,*) nseed
4828 write(icsa_history,*) ntotal,cut1,cut2
4829 write(icsa_history,*) estop
4830 write(icsa_history,*) icmax,irestart
4831 write(icsa_history,*) ntbankm,dele,difcut
4832 write(icsa_history,*) iref,rmscut,pnccut
4833 write(icsa_history,*) ndiff
4835 write(icsa_history,*)
4838 open(icsa_bank1,file=csa_bank1,status="unknown")
4839 write(icsa_bank1,*) 0
4843 end subroutine initial_write
4844 !-----------------------------------------------------------------------------
4845 subroutine restart_write
4848 ! implicit real*8 (a-h,o-z)
4849 ! include 'DIMENSIONS'
4850 ! include 'COMMON.IOUNITS'
4851 ! include 'COMMON.CSA'
4852 ! include 'COMMON.BANK'
4854 ! integer :: ntf,ik,iw_pdb
4855 ! real(kind=8) :: var,ene
4857 #if defined(AIX) || defined(PGI)
4858 open(icsa_history,file=csa_history,position="append")
4860 open(icsa_history,file=csa_history,access="append")
4862 write(icsa_history,*)
4863 write(icsa_history,*) "This is restart"
4864 write(icsa_history,*)
4865 write(icsa_history,*) nconf
4866 write(icsa_history,*) jstart,jend
4867 write(icsa_history,*) nstmax
4868 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4869 write(icsa_history,*) nran0,nran1,irr
4870 write(icsa_history,*) nseed
4871 write(icsa_history,*) ntotal,cut1,cut2
4872 write(icsa_history,*) estop
4873 write(icsa_history,*) icmax,irestart
4874 write(icsa_history,*) ntbankm,dele,difcut
4875 write(icsa_history,*) iref,rmscut,pnccut
4876 write(icsa_history,*) ndiff
4877 write(icsa_history,*)
4878 write(icsa_history,*) "irestart is: ", irestart
4880 write(icsa_history,*)
4884 end subroutine restart_write
4885 !-----------------------------------------------------------------------------
4887 !-----------------------------------------------------------------------------
4888 subroutine write_pdb(npdb,titelloc,ee)
4890 ! implicit real*8 (a-h,o-z)
4891 ! include 'DIMENSIONS'
4892 ! include 'COMMON.IOUNITS'
4893 character(len=50) :: titelloc1
4894 character*(*) :: titelloc
4895 character(len=3) :: zahl
4896 character(len=5) :: liczba5
4898 integer :: npdb !,ilen
4902 ! real(kind=8) :: var,ene
4906 if (npdb.lt.1000) then
4907 call numstr(npdb,zahl)
4908 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
4910 if (npdb.lt.10000) then
4911 write(liczba5,'(i1,i4)') 0,npdb
4913 write(liczba5,'(i5)') npdb
4915 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
4917 call pdbout(ee,titelloc1,ipdb)
4920 end subroutine write_pdb
4921 !-----------------------------------------------------------------------------
4923 !-----------------------------------------------------------------------------
4924 subroutine write_thread_summary
4925 ! Thread the sequence through a database of known structures
4926 use control_data, only: refstr
4928 use energy_data, only: n_ene_comp
4930 ! implicit real*8 (a-h,o-z)
4931 ! include 'DIMENSIONS'
4933 use MPI_data !include 'COMMON.INFO'
4936 ! include 'COMMON.CONTROL'
4937 ! include 'COMMON.CHAIN'
4938 ! include 'COMMON.DBASE'
4939 ! include 'COMMON.INTERACT'
4940 ! include 'COMMON.VAR'
4941 ! include 'COMMON.THREAD'
4942 ! include 'COMMON.FFIELD'
4943 ! include 'COMMON.SBRIDGE'
4944 ! include 'COMMON.HEADER'
4945 ! include 'COMMON.NAMES'
4946 ! include 'COMMON.IOUNITS'
4947 ! include 'COMMON.TIME1'
4949 integer,dimension(maxthread) :: ip
4950 real(kind=8),dimension(0:n_ene) :: energia
4952 integer :: i,j,ii,jj,ipj,ik,kk,ist
4953 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
4955 write (iout,'(30x,a/)') &
4956 ' *********** Summary threading statistics ************'
4957 write (iout,'(a)') 'Initial energies:'
4958 write (iout,'(a4,2x,a12,14a14,3a8)') &
4959 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
4960 'RMSnat','NatCONT','NNCONT','RMS'
4961 ! Energy sort patterns
4966 enet=ener(n_ene-1,ip(i))
4969 if (ener(n_ene-1,ip(j)).lt.enet) then
4971 enet=ener(n_ene-1,ip(j))
4983 ist=nres_base(2,ii)+ipatt(2,i)
4985 energia(i)=ener0(kk,i)
4987 etot=ener0(n_ene_comp+1,i)
4988 rmsnat=ener0(n_ene_comp+2,i)
4989 rms=ener0(n_ene_comp+3,i)
4990 frac=ener0(n_ene_comp+4,i)
4991 frac_nn=ener0(n_ene_comp+5,i)
4994 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4995 i,str_nam(ii),ist+1,&
4996 (energia(print_order(kk)),kk=1,nprint_ene),&
4997 etot,rmsnat,frac,frac_nn,rms
4999 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5000 i,str_nam(ii),ist+1,&
5001 (energia(print_order(kk)),kk=1,nprint_ene),etot
5004 write (iout,'(//a)') 'Final energies:'
5005 write (iout,'(a4,2x,a12,17a14,3a8)') &
5006 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5007 'RMSnat','NatCONT','NNCONT','RMS'
5011 ist=nres_base(2,ii)+ipatt(2,i)
5013 energia(kk)=ener(kk,ik)
5015 etot=ener(n_ene_comp+1,i)
5016 rmsnat=ener(n_ene_comp+2,i)
5017 rms=ener(n_ene_comp+3,i)
5018 frac=ener(n_ene_comp+4,i)
5019 frac_nn=ener(n_ene_comp+5,i)
5020 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5021 i,str_nam(ii),ist+1,&
5022 (energia(print_order(kk)),kk=1,nprint_ene),&
5023 etot,rmsnat,frac,frac_nn,rms
5025 write (iout,'(/a/)') 'IEXAM array:'
5026 write (iout,'(i5)') nexcl
5028 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5030 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5031 'Max. time for threading step ',max_time_for_thread,&
5032 'Average time for threading step: ',ave_time_for_thread
5034 end subroutine write_thread_summary
5035 !-----------------------------------------------------------------------------
5036 subroutine write_stat_thread(ithread,ipattern,ist)
5038 use energy_data, only: n_ene_comp
5040 ! implicit real*8 (a-h,o-z)
5041 ! include "DIMENSIONS"
5042 ! include "COMMON.CONTROL"
5043 ! include "COMMON.IOUNITS"
5044 ! include "COMMON.THREAD"
5045 ! include "COMMON.FFIELD"
5046 ! include "COMMON.DBASE"
5047 ! include "COMMON.NAMES"
5048 real(kind=8),dimension(0:n_ene) :: energia
5050 integer :: ithread,ipattern,ist,i
5051 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
5053 #if defined(AIX) || defined(PGI)
5054 open(istat,file=statname,position='append')
5056 open(istat,file=statname,access='append')
5059 energia(i)=ener(i,ithread)
5061 etot=ener(n_ene_comp+1,ithread)
5062 rmsnat=ener(n_ene_comp+2,ithread)
5063 rms=ener(n_ene_comp+3,ithread)
5064 frac=ener(n_ene_comp+4,ithread)
5065 frac_nn=ener(n_ene_comp+5,ithread)
5066 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5067 ithread,str_nam(ipattern),ist+1,&
5068 (energia(print_order(i)),i=1,nprint_ene),&
5069 etot,rmsnat,frac,frac_nn,rms
5072 end subroutine write_stat_thread
5073 !-----------------------------------------------------------------------------
5075 !-----------------------------------------------------------------------------
5076 end module io_config